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Abstract 

Diode  lasers  are  useful  in  military  and  commercial  applications  that  have  strict 
requirements  for  size,  weight  and  power.  This  includes  the  use  of  diode  lasers  in 
optoelectronic  and  photonic  integrated  circuits,  which  can  lead  to  new  technologies 
in  optical  communications  and  optical  interconnects  in  high  performance  computing 
systems.  For  these  systems  to  be  effective,  the  diode  laser  must  be  modulated  at 
frequencies  beyond  current  limits  which  are  typically  a  few  GHz.  This  barrier  can  be 
broken  by  optically  coupling  a  diode  laser  with  a  similar  laser. 

A  set  of  single  mode  rate  equations  models  the  dynamics  of  twin  optically  cou¬ 
pled  diode  lasers.  Steady  states  of  the  system  are  derived  analytically  or  calculated 
numerically  when  an  analytic  expression  is  not  easily  available.  The  stability  of  the 
steady  states  is  examined  by  using  a  linear  stability  analysis,  which  is  also  used  in 
an  algorithm  that  calculates  the  infinitesimal  modulation  response.  The  modulation 
response  is  also  calculated  by  using  a  numerical  method  that  directly  integrates  the 
rate  equations.  Typical  parameters  for  an  InGaAsP  diode  laser  are  used  in  the  algo¬ 
rithms  to  investigate  mutual  coupling  and  evanescent  coupling.  It  will  be  shown  that 
mutually  coupled  lasers  can  be  effectively  modulated  out  to  frequencies  of  approxi¬ 
mately  9  GHz  compared  to  4  GHz  for  a  solitary  laser.  For  evanescent  coupling,  the 
steady  states  are  unstable  over  large  regions  of  the  parameter  space,  but  this  is  reme¬ 
died  by  introducing  asymmetric  DC  currents  through  the  lasers  or  by  introducing  the 
effects  of  gain  saturation.  With  stable  steady  states,  evanescently  coupled  lasers  can 
be  effectively  modulated  at  frequencies  out  to  about  30  GHz  which  is  more  than  a 
seven  fold  improvement  over  a  solitary  laser. 
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MODULATION  RESPONSE 


OF 

TWIN  OPTICALLY  COUPLED 
DIODE  LASERS 

I.  Introduction 

Diode  lasers  are  useful  in  high  speed  optical  communications,  as  well  as  long-haul 
terrestrial  and  transoceanic  communications  [1;  2],  Because  of  their  small  size  and 
high  efficiency,  they  are  also  useful  in  UAV  and  satellite  payloads,  as  well  as  in  other 
remote  tactical  applications  that  have  strict  requirements  on  size,  weight  and  power 
[1].  Diode  lasers  undergoing  optical  feedback  or  injection  can  be  driven  into  chaos. 
In  this  state  diode  lasers  can  be  used  for  secure  transmission  of  information  [3;  4], 
Diode  lasers  can  also  play  an  important  role  in  optoelectronic  and  photonic  integrated 
circuits  as  small  efficient  sources  of  photons  [5].  These  optoelectronic  and  photonic 
devices  can  possibly  lead  to  the  next  generation  of  communications  and  computing 
devices,  where  the  circuitry  makes  use  of  photons  instead  of  electrons  to  process 
information  and  perform  computations. 

In  order  to  transmit  information,  the  output  of  the  diode  laser  must  be  mod¬ 
ulated.  This  is  commonly  achieved  in  diode  lasers  by  varying  the  electrical  current 
through  the  laser  [6].  The  diode  laser  must  be  effectively  modulated  at  high  fre¬ 
quencies  for  fast  transmission  of  information.  While  diode  lasers  are  small,  light  and 
efficient,  they  typically  do  not  perform  effectively  when  modulated  at  frequencies 
beyond  the  relaxation  oscillation  frequency,  which  is  typically  a  few  gigahertz  [6]. 
However,  it  has  been  shown  in  experiments  [7;  8]  and  theory  [1;  2;  6]  that  diode  lasers 
can  be  effectively  modulated  beyond  the  frequency  of  relaxation  oscillations  when 
optically  injected  or  optically  coupled  with  similar  diode  lasers.  Because  of  this,  the 
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modulation  characteristics  of  coupled  diode  lasers  have  become  a  topic  of  interest 
[7;  8;  9;  10;  11;  12]. 


The  modulation  characteristics  of  a  diode  laser  are  examined  by  studying  the 
modulation  response.  The  modulation  response  of  a  diode  laser  is  defined  as  the 
ratio  of  the  modulated  electric  held  amplitude  to  the  modulated  electrical  current 
amplitude  through  the  laser  normalized  to  the  response  at  a  modulation  frequency 
of  zero  [1;  13].  From  the  modulation  response,  characteristics  such  as  bandwidth 
and  peak  response  are  extracted.  Figure  1  shows  some  features  of  the  modulation 
response  of  a  solitary  diode  laser.  For  communications  and  other  applications,  the 
bandwidth  should  be  broad,  which  allows  more  information  to  be  transmitted,  and 
the  modulation  response  as  hat  as  possible  across  the  bandwidth  in  order  to  avoid 
having  the  band  dominated  by  one  or  possibly  a  few  modulation  frequencies. 


Figure  1.  Some  properties  of  the  modulation  response  for  a  solitary  InGaAsP  diode 
laser.  Adapted  from  Reference  [1]  and  Reference  [14]. 
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In  order  to  calculate  the  modulation  response,  the  dynamics  of  the  laser  must 
be  studied.  The  first  step  is  to  derive  [13;  14;  15]  an  appropriate  set  of  rate  equations. 
For  a  solitary  laser,  a  set  of  equations  that  describe  the  loss  and  gain  of  photons  as 
well  as  the  loss  and  gain  of  electrons  in  the  gain  medium  is  sufficient  [14].  However, 
the  model  of  a  coupled  laser  system  must  account  for  the  electric  field  amplitude  and 
phase,  since  the  electric  fields  of  the  twin  lasers  must  be  coherently  added  [16].  In 
Section  2.1,  a  rate  equation  for  the  complex  slowly  varying  electric  held  is  formulated 
and  used  in  Section  2.2  to  develop  a  model  for  a  system  of  twin  optically  coupled  diode 
lasers  [11;  12].  The  model  is  then  applied  to  a  solitary  laser  in  Section  2.3.  This  will 
establish  a  performance  baseline  for  the  coupled  system  as  well  as  define  some  terms 
commonly  used  in  this  held  of  study.  The  treatment  used  for  the  solitary  laser  is 
then  used  in  Chapter  III  to  develop  a  set  of  computational  methods  that  determine 
the  stability  of  the  steady  states  (Section  3.1),  compute  the  inhnitesimal  modulation 
response  (Section  3.2),  and  numerically  compute  the  modulation  response  (Section 
3.3).  Next,  these  computational  methods  are  used  in  Chapter  IV  to  simulate  two 
cases  of  optical  coupling  (Figure  2).  First  in  Section  4.1,  where  the  lasers  are  mutually 
coupled  (Figure  2a)  by  directly  injecting  photons  through  the  output  couplers  [17] 
and  second  where  photons  are  passed  through  the  sides  of  the  lasers  (Figure  2b)  by 
way  of  an  evanescent  wave  [1;  6;  18;  19].  It  is  shown  in  this  document  that  the 
bandwidth  of  mutually  coupled  lasers  is  slightly  enhanced  when  modulated  in-phase 
with  the  same  signal,  but  when  modulated  out-of-phase  the  response  becomes  flatter. 
It  is  also  shown  in  this  document  that  the  modulation  bandwidth  of  evanescently 
coupled  diode  lasers  is  about  six  times  that  of  a  solitary  laser  when  the  coupled  lasers 
are  modulated  out-of-phase  with  the  same  signal. 
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(a)  Mutually  Coupled  Diode  Lasers 


(b)  Evanescently  Coupled  Diode  Lasers 


Laser  1  Active  Region 

(2)  Laser  2  Active  Region 

(3)  Substrate 

(4)  Metal  Contact 

(5)  N-Type  Confinement  Layer 

(6)  P-Type  Confinement  Layer 


Figure  2.  (a)  Mutually  coupled  buried-heterostructure  twin  diode  lasers.  Between 

the  lasers  the  photons  couple  out  of  the  output  couplers  and  into  the  output  couplers 
of  the  other  laser,  (b)  Evanescently  coupled  buried-heterostructure  twin  diode  lasers. 
The  confinement  layers  set  up  the  lasing  mode  through  total  internal  reflection.  When 
the  lasers  are  close  enough  together  the  evanescent  waves  in  each  laser  will  couple 
into  the  lasing  mode  of  the  other.  Adapted  from  Reference  [13]  and  Reference  [20]. 
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II.  Theory 

In  this  document  single  mode  rate  equations  are  used  to  model  the  dynamics  of  twin 
optically  coupled  diode  lasers.  Because  of  this,  it  is  important  to  understand  how 
these  rate  equations  are  derived.  The  first  step  is  to  examine  a  standard  set  of  rate 
equations  for  solitary  diode  lasers  [14;  21].  These  standard  equations  describe  the 
photon  and  carrier  densities  in  the  diode  laser,  but  in  order  to  model  an  optically 
coupled  system,  rate  equations  for  the  electric  field  amplitude  and  phase  must  be 
derived  [13].  Once  this  is  complete,  the  model  for  a  coupled  system  is  developed. 
Finally,  the  model  is  applied  to  a  solitary  laser,  which  establishes  a  performance 
baseline  and  leads  into  the  development  of  computational  methods  for  the  coupled 
system. 

2. 1  Rate  Equations 

Diode  laser  systems  are  modeled  by  a  set  of  rate  equations  that  describe  the 
exchange  of  carriers  and  photons  [14;  21].  The  rate  equation  for  the  photon  density 

is 


—  =  G(N,S)S(t)  +  Rap-^-7  (1) 

Ctt  Tp 

and  includes  terms  that  account  for  gains  in  photons  due  to  spontaneous  and  stimu¬ 
lated  emission  and  losses  due  to  absorption  and  output  coupling.  The  rate  equation 
for  the  carrier  density  is 

^-  =  P(t)-^-G(N,S)S{t),  (2) 

at  tc 

and  includes  terms  that  account  for  gain  due  to  pumping  and  losses  due  to  stimulated 
emission,  spontaneous  emission  and  non-radiative  transitions.  In  Equations  (1,2)  t 
is  time  in  seconds,  S  is  the  photon  density  in  cm-3,  G  is  net  gain  per  second,  and 
Rsp  is  the  rate  of  photons  that  are  spontaneously  emitted  into  the  lasing  mode  in 
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cnr3s_1.  The  term  that  includes  the  photon  lifetime  in  seconds  (rp)  accounts  for 
the  loss  of  photons  due  to  absorption  and  output  coupling.  The  term  that  includes 
the  carrier  density  in  cm”3  (IV)  accounts  for  the  number  of  electron-hole  pairs  in 
the  active  region.  The  electrical  pumping  rate  in  cm”3s_1  (P)  is  proportional  to  the 
current  density  through  the  diode  laser.  The  term  that  includes  the  carrier  lifetime 
in  seconds  (rc)  accounts  for  electron-hole  recombination  for  both  radiative  and  non- 
radiative  transitions.  Figure  3  illustrates  the  processes  in  a  diode  laser  that  are 
described  by  the  rate  equations. 


Pit) 

Conduction  *  # 

Band  • 


Figure  3.  Lasing  in  a  diode  laser.  Gray  arrows  indicate  photons  that  contribute 
to  the  lasing  mode.  Solid  black  circles  indicate  electrons  (carriers).  Hollow  black 
circles  indicate  holes.  The  large  gray  circle  indicates  spontaneously  emitted  photons 
in  the  lasing  mode.  The  large  hollow  circle  indicates  spontaneously  emitted  photons 
not  in  the  lasing  mode  and  the  curved  black  line  indicates  non-radiative  transitions. 
Downward  pointing  black  arrows  indicate  electron  hole  recombination  and  the  upward 
pointing  black  arrow  indicates  electron  promotion  to  the  conduction  band  through 
photon  absorption. 


If  the  photon  density  is  assumed  to  be  at  a  level  where  gain  saturation  effects 
are  minimal,  then  a  linear  model  of  gain  [14;  22]  can  be  used.  The  linear  gain  model 
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IS 


G(N)  =  Gth  +  Gn(N  -  Nth) 


1 


+  Gn(N  —  Nth), 


(3) 


where  Gth  is  the  threshold  gain  with  units  of  s  1,  G^  is  the  differential  gain  with 


units  of  cm3/s,  and  Nth  is  the  threshold  carrier  density.  Since  accounts  for  all 


losses  in  the  lasing  mode,  the  threshold  gain  equals  W  If  the  number  of  photons  in 
the  lasing  mode  is  primarily  due  to  stimulated  emission,  which  is  typically  the  case 
while  operating  above  threshold,  then  the  rate  of  spontaneous  emission  into  the  lasing 
mode  can  be  ignored.  With  this  assumption  and  the  linear  gain  model,  Equations 
(1,2)  become 


—  =  GN(N-Nth)S{t ) 


(4) 


(5) 


This  set  of  rate  equations  can  be  used  to  model  the  dynamics  of  a  solitary  diode  laser, 
but  in  order  to  model  coupled  diode  lasers  the  fields  created  by  stimulated  emission 


and  the  fields  created  by  coupling  within  each  laser  must  be  coherently  added.  This 


requires  rate  equations  for  the  electric  fields  in  each  laser.  The  place  to  start  for  these 
rate  equations  is  Maxwell’s  wave  equation  [13;  14;  15]. 

The  following  derivation  follows  Agrawal  and  Dutta’s  treatment  in  Reference 
[13].  Assuming  propagation  along  an  arbitrary  z-axis,  non-conductive  and  nonmag¬ 
netic  materials,  instantaneous  changes  in  the  polarization  held,  and  single  mode  op¬ 
eration,  the  wave  equation  for  the  electric  held  is  written  as 


(6) 
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where  £  is  the  electric  field  in  Gaussian  units,  c  is  the  speed  of  light  in  a  vacuum, 
and  e  is  the  relative  permittivity.  The  electric  field  is  written  as 

£ (x,  y,  z,  t )  =  ip(x,  y )  sin (kz)E(t)e~luJt,  (7) 


where  E  is  the  slowly  varying  electric  field,  u  is  the  angular  frequency  of  the  lasing 
mode,  and  k  is  the  effective  mode  propagation  constant.  The  transverse  electric 
field  distribution  is  represented  by  ^>{x,y)  while  sin (kz)  accounts  for  the  electric  field 
distribution  along  the  direction  of  propagation  created  by  the  Fabry-Perot  cavity 
[13].  Substituting  Equation  (7)  into  Equation  (6)  and  ignoring  second  derivatives  of 
E,  since  it  is  assumed  to  be  slowly  varying,  yields 

(n  ■  ap  2  \ 

— ~dt  ~  'E^J  ' 

The  spatial  dependence  of  Equation  (8)  is  handled  by  multiplying  Equation  (8)  by 
ip(x,y)  sin(kz)  and  integrating  across  all  x,  y,  and  z.  Assuming  that  ip(x,  y)  sin(fcz) 
is  properly  normalized  such  that 


^(x,  y)2  s'm2 (kz)dxdydz  =  1, 


(9) 


Equation  (8)  becomes 


<e> 


2 iu  dE 
c2  dt 


LU 


=  (c)- o-fc2  m 


where 


(10) 


<£>  = 


(x,  y)^(x:  y)2  sin 2{kz)dxdydz. 


Since  the  electric  field  will  most  likely  be  dispersed  across  several  materials,  e  may  not 
be  constant  across  the  laser.  This  is  the  reason  for  (e),  which  is  a  weighted  average 
of  relative  permittivities  across  the  areas  where  there  is  an  electric  field  and  can  be 


approximated  as  the  square  of  the  effective  mode  index.  While  lasing,  the  carrier 
density  is  pinned  at  its  threshold  value  as  well  as  the  gain.  This  means  that  the 
relative  permittivity  is  real  at  threshold,  and  any  deviations  in  the  carrier  density 
will  lead  to  a  change  in  relative  permittivity  [15].  This  is  accounted  for  with 

(e)  =  (nr  +  Snr  +  iSrii)2 
~  n*  +  2inr8rii(l  —  ia), 


where 


8nr 


(12) 


In  Equation  (11),  nr  is  the  real  part  of  the  effective  mode  index,  n,  is  the  imaginary 
part  of  the  effective  mode  index,  8  indicates  small  deviations,  and  a  is  the  line-width 
enhancement  factor  [15].  Since 


k 


(13) 


Equation  (10)  is  written  as 


d  E 
~dt 


u)8rii 

nr 


E(t). 


(14) 


The  imaginary  part  of  the  refractive  index  is  related  to  net  gain  through 

=  ,16) 

With  Equation  (15)  and  using  the  linear  gain  model  in  Equation  (3),  Equation  (14) 
is  written  as 


nr  i 

—  =  -(l  -ia)GN(N-Nth)E(t), 


(16) 
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which  is  the  rate  equation  for  the  slowly  varying  electric  field.  Since  the  square  of  the 
absolute  value  of  the  slowly  varying  electric  field  is  proportional  to  the  intensity  and 
intensity  is  proportional  to  the  photon  density,  the  electric  field  can  be  normalized 
such  that  S  =  EE*.  With  this  normalization,  the  rate  equations  for  the  slowly 
varying  electric  field  and  carrier  density  become 

^  =  1(1  _  ia)GN{N  -  Nth)E(t)  (17) 

at  2 

cllL  =  P[t)  _  M  _  (I  +  Gn(N  -  Nth))\E(t)\ 2.  (18) 

Cit  Tc  Tp 

Table  1  shows  some  typical  values  of  the  parameters  in  Equations  (17,  18).  These 
equations  can  be  used  to  model  a  solitary  laser,  but  an  additional  term  will  need  to 
be  added  to  Equation  (17)  in  order  to  model  diode  lasers  undergoing  optical  injection 
or  feedback.  This  is  converted  in  the  next  section  where  the  model  for  twin  optically 
coupled  diode  lasers  is  developed. 


Table  1.  Parameters  for  a  1.3  pm  buried-heterostructure  InGaAsP  diode  laser  [13]. 


Name 

Symbol  Value 

Units 

Line-width  Enhancement  Factor 

a 

5 

- 

Carrier  Lifetime 

Tc 

2.2 

ns 

Photon  Lifetime 

Tp 

1.6 

ps 

Time  in  Cavity 

1 in 

2.8 

ps 

Differential  Gain 

Gn 

6.6  x  10~7 

cm3/s 

Threshold  Carrier  Density 

Nth 

1.0  x  1018 

cm-3 

Threshold  Current 

- 

15.8 

mA 

Cavity  Length 

- 

250 

/irn 

Active  Region  Width 

- 

2 

pm 

Active  Region  Thickness 

- 

0.2 

pm 

Effective  Mode  Index 

- 

3.4 

- 

2.2  Model  of  Twin  Optically  Coupled  Diode  Lasers 

A  diode  laser  undergoing  optical  injection  or  feedback  is  modeled  by  a  modified 
version  of  the  rate  equations  derived  in  Section  2.1.  Lang  and  Kobayashi  show  that 
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the  rate  equation  for  the  electric  field  needs  an  additional  term  in  order  to  account 
for  optical  injection  or  feedback  [16].  The  modified  rate  equations  are 


=  1  -  ioi)GN(N  -  Nth)E(t)  +  —Einj(t  -  t )  (19) 

dt  Z  7~in 

E  =  P[t)  -  M  -  (I  +  Gn(N  -  Nih))\E(t)\\  (20) 

Cit  Tc  Tp 


where  k  is  the  optical  coupling  term,  T;n  is  the  time  that  a  photon  takes  to  transverse 
the  laser  cavity  twice,  and  r  is  the  time  delay  of  the  optical  injection  or  feedback. 
For  convenience,  Equations(19,20)  are  transformed  into  dimensionless  form  [22;  23] 
by  introducing 


E  =  sJ^^E 
N  =  ^(N-Nth) 


p  —  EpG^Nth  f  P  _ 


P  ,  —  Nth 
*  th  _ 


p  —  Tc 
~  TP 
~  _  KTV 

K  =  — - 


V  Eh 


T  =  — 


(21) 


In  these  equations,  E  is  the  dimensionless  slowly  varying  electric  held,  N  is  the 
dimensionless  excess  carrier  density,  t  is  dimensionless  time,  P  is  the  dimensionless 
excess  electrical  pumping  rate,  Pth  is  the  threshold  electrical  pumping  rate,  and  T  is 
the  ratio  of  the  carrier  lifetime  to  the  photon  lifetime.  Table  2  shows  some  values  of 
the  coefficients  in  Equation  (21)  using  values  from  Table  1.  The  rate  equations  in 
terms  of  the  new  variables  in  Equation  (21)  now  become 


dE 

dt 

dN 

dt 


(1  -  ia)N(i)E(t )  +  kEinj{t  -  t) 
P(t)  -  N(i)  -  (1  +  2N(t))\E(t)\\ 


(22) 

(23) 
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Table  2.  Dimensionless  conversion  coefficients  using  values  from  Table  1. 


Coefficient 

Value 

Units 

\J  tcGn/2 

2.7  x  10"8 

cm3/2 

tpGn/^ 

5.3  x  KT19 

cm3 

TpGNNth/2 

0.53 

— 

T 

1375 

— 

7~p/ En 

0.57 

— 

For  a  twin  optically  coupled  diode  laser  system  each  laser  is  modeled  by  Equa¬ 
tions  (22,23)  with  the  output  of  each  laser  injecting  into  the  other  as  shown  in  the 
following  equations: 


dE 

— =  (1  —  ia)NiE\  +  kE2(t  —  t)  (24) 

dt 

El  =  (1  -  ia)N2E2  +  kEi(t  -  t)  (25) 

dt 

tE±  =  p1  -Ni  -  (l  +  2iV1)|E1|2  (26) 

dt 

TEl  =  P2-N2-(1  +  2N2)\E2\2.  (27) 

dt 


Since  both  E  and  k  are  complex,  Equations  (24-27)  are  decomposed  into  real  and 
imaginary  parts  by  substituting 


E(t)  =  R(i)ea® 
k  —  r}e 


(28) 

(29) 


where  R  is  the  dimensionless  slowly  varying  electric  held  amplitude,  6  is  the  slowly 
varying  electric  held  phase,  77  is  the  dimensionless  coupling  strength,  and  (j)  is  the 
coupling  phase.  After  substitution  Equations  (24,25)  are  decomposed  into  their  am- 
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plitude  and  phase  parts  as  shown  in  the  following  equations: 


,JD 

— J-  =  N1R1  +  r]R2(i  -  t)  cos (d2(t  -  t)  -  91  +  0)  (30) 

at 

=  N2R2  +  r]Ri(t  -  t)  cos(6>i (t  -  t )  -  02  +  0)  (31) 

at 

=  -a>N i  +r]R2^  T ^  sin (02(t  -  t)  -  6 1  +  0)  (32) 

at 

=  -aN2  +  — —  sin(6*i(t  —  t)  —  02  +  (j>).  (33) 

dt  ht  2 


This  model  is  simplihed  by  only  considering  instantaneous  coupling  where  f  =  0.  The 
slowly  varying  electric  held  phase  difference, 

e  =  02-0u  (34) 

is  now  introduced  to  reduce  the  system  of  six  equations  to  a  system  of  five  equations 
as  shown  in  the  following  equations: 

,  7  r> 

— J-  =  N±Ri  +  ijR2  cos(©  +  0)  (35) 

dt 

— ^  =  N2R2  +  r)R,]  cos(©  —  0)  (36) 

dt 

_  =  —a(N2  —  Ni)  —r)  ( sin(©  -  0)  +  sin(©  +  0)  J  (37) 

rdiVi  =  ^  ^  +  2Nl)R\  (38) 

dt 

T-jf-  =  A  ~  N2  -  (1  +  2 N2)Rl  (39) 

dt 

Table  3  shows  some  approximate  values  of  the  parameters  in  Equations  (35-39).  This 
system  of  equations  along  with  the  appropriate  computational  methods  are  used  to 
investigate  the  modulation  response  of  evanescently  coupled  twin  diode  lasers  and 
mutually  coupled  twin  diode  lasers,  but  an  investigation  of  the  modulation  response 
of  a  solitary  laser  is  provided  in  the  next  section  to  establish  a  performance  baseline 
and  to  define  some  terms. 
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Table  3.  Typical  values  of  the  parameters  in  Equations  (35-39). 


Symbol 

Typical  Value 

Reference 

a 

Rj  5 

(1]  [6]  [15]  [22] 

rj 

Rj  10"2 

[22] 

T 

Rj  103 

[1][6]  [22] 

P 

Rj  1 

[1]  [6]  [22] 

2.3  Modulation  Response  of  a  Solitary  Diode  Laser 

Following  the  procedure  found  in  Reference  [14]  and  [21],  The  rate  equations 
derived  in  Section  2.2  are  used  to  model  a  solitary  diode  laser  and  solve  for  its  in¬ 
finitesimal  modulation  response.  This  lays  out  a  performance  baseline  for  the  twin 
optically  coupled  diode  laser  system,  defines  some  terms  and  outlines  the  basic  pro¬ 
cedure  for  the  linear  stability  analysis  that  will  be  used  to  analyze  the  coupled  laser 
system. 

The  solitary  laser  is  modeled  by  taking  rj  —  0  in  Equations  (35-39)  and  choosing 
one  set  of  subscripts  to  arrive  at 


dR 

dt 

dd_ 

dt 

TdN 

dt 


NR 

-aN 

P  -  N  -  (1  +  2N)R2. 


(40) 

(41) 

(42) 


To  End  the  infinitesimal  modulation  response  Equations  (40-42)  are  linearized  about 
their  steady  state  values.  The  steady  state  values  are  found  by  equating  ^  in  Equa¬ 
tion  (40)  to  zero.  This  means  that  either  R  =  0  or  N  =  0.  Since  R  p-  0,  otherwise 
the  laser  is  off,  then  N  —  0,  which  means  the  laser  is  operating  at  threshold.  This 
result  is  then  substituted  into  Equation  (41)  and  Equation  (42)  to  find  the  steady 
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state  values  of  R  and  9.  The  steady  state  solutions  are 


Ns  =  0  (43) 

Rs  =  \fps ,  (44) 

where  Ns  is  the  steady  state  dimensionless  excess  carrier  density,  Rs  is  the  steady 
state  dimensionless  slowly  varying  electric  field  amplitude,  and  Ps  is  the  steady  state 
dimensionless  excess  electrical  pumping  rate.  The  steady  state  slowly  varying  electric 
field  phase  (9S)  becomes  an  arbitrary  constant.  The  linearization  is  accomplished  by 
substituting 


R(t)  =  RS  +  6R(t)  (45) 

9(t)  =  9S  +  69(t)  (46) 

N{t)  =  NS  +  6N(t )  (47) 

P(t)  =  PS  +  6P(t)  (48) 


into  Equations  (40-42).  The  equations  are  then  expanded  while  eliminating  terms 
such  as  5R2  and  5R5N  since  the  6  terms  are  defined  as  small  oscillations  about  the 
steady  state  values.  The  results  of  this  linearization  are 


l.SR 

dt 

ise 

dt 

T—JN 

dt 


—aSN 

-(1  +  2 PS)5N  -  2  JPsSR  +  5P. 


(49) 

(50) 

(51) 


This  set  of  equations  can  also  be  written  in  matrix  form, 


SR 

0 

0 

\fPs 

6R 

0 

69 

= 

0 

0 

—a 

69 

+ 

0 

6N 

-2  \[PsjT 

0 

-(1  +  2 Ps)/T  _ 

6N 

1/T 

SP, 


(52) 
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where  the  3x3  matrix  is  the  Jacobian  of  Equations  (40-42)  when  the  dimensionless 
excess  electrical  pumping  rate  is  held  constant.  The  dynamics  of  this  system  will  now 
be  examined  with  no  modulation  signal,  5P  =  0,  in  order  to  define  some  terms  that 
will  be  used  throughout  this  document. 

Taking  5P  =  0,  Differentiating  Equation  (51)  and  substituting  Equation  (49) 
for  gives  a  second  order  ordinary  differential  equation. 

„  1  _i_  2P  ,  9 P 

5N"  +  T  s  5 N'  +  -j^5N  =  0.  (53) 

The  solution  of  Equation  (53)  is 

SN{t)  =  Ciex+i  +  C2ex~\  (54) 


where 


A+ 

A_ 


1  +  2  Ps) 

+  i\ 

2  Ps 

2  T  J 

T 

1  +  2  Ps\ 

-i. 

2  Ps 

2  T  J 

T 

1  +  2  Ps\ 

2  T  J 


1  +  2  Ps\ 
2  T  J 


are  the  eigenvalues  of  the  Jacobian  matrix  in  Equation  (52)  and  Cn  are  arbitrary 
constants.  This  solution  is  a  damped  simple  harmonic  oscillator  with  a  dimensionless 
relaxation  damping  coefficient, 


_  1  +  2PS 

~/r  -  2T 


(55) 


and  a  dimensionless  relaxation  oscillation  frequency, 


(56) 
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This  oscillation  frequency  is  approximated  as 


(57) 

since  T  is  approximately  103  making  2 Ps/T  »  7)?.  The  approximate  form  of  the 
dimensionless  relaxation  oscillation  frequency  in  Equation  (57)  will  be  used  when 
referring  to  this  frequency.  These  terms  can  now  be  used  in  the  derivation  of  the  in- 
finitessimal  modulation  response  as  well  as  in  subsequent  equations  in  this  document. 

The  infinitesimal  modulation  response  is  found  by  taking  the  Fourier  transform 
of  Equations  (49,53)  as  shown  in  the  following  equations: 


=  \[^SR{5N} 

=  -u>lF{S  N}  +  2  lrujmT{5N)  +  w2rT{6N}, 


(58) 

(59) 


where  JF{}  indicates  a  Fourier  transform  operation  and  a ;m  is  the  dimensionless  angu¬ 
lar  modulation  frequency.  Equations  (58,59)  are  now  solved  for  the  Fourier  transform 
of  the  electric  field  divided  by  the  Fourier  transform  of  the  carrier  density  to  arrive 
at 


T{SR]  _  ^f¥s 

R{5P}  T(uj 2  -  a 4  +  2i^rujm) ' 

Taking  the  absolute  value  yields 


(60) 


R{5R} 

R{5P} 


VPs _ 

-wr2) 2  +  47 2ruRm 


(61) 


Finally,  Equation  (61)  is  normalized  to  the  response  at  zero  frequency  to  arrive  at 


M(ujm ) 


u 


2 

r 


Wr)2  +  47r2a;m 


where  M  denotes  the  modulation  response. 


(62) 
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The  infinitesimal  modulation  response  is  also  found  by  using  the  matrix  form  of 
the  linearized  equations.  This  treatment  is  the  basis  for  calculating  the  infinitesimal 
modulation  response  of  the  coupled  system.  The  first  step  is  taking  the  Fourier 
transform  of  Equation  (52)  to  arrive  at 


R{8R} 

R{8R} 

0 

iiOm 

T{89} 

F{8N} 

=  ,J 

T{59} 

T{8N} 

+ 

0 

1/T 

R{8P}, 

(63) 

where  J  is  the  Jacobian  matrix.  Collecting  like  terms  and  dividing  by  the  Fourier 
transform  of  the  electrical  pumping  rate  modulation  yields 


•FjM} 

0 

Hs  p} 

CF{59} 

Hsp) 

[t / 

0 

P{SN} 

1 

L  T{5P}  J 

L  T  J 

(64) 


where  /  is  the  identity  matrix.  Finally,  the  modulation  response  is  found  by  multi¬ 
plying  the  first  row  third  column  of  the  matrix  in  Equation  (64)  by  —7^,  taking  the 
absolute  value,  and  normalizing  to  the  response  at  a  modulation  frequency  of  zero. 


[J 


ICO, 


J] 


-1 

1,3 


Ji 


-1 


(65) 


Since  J  is  known  analytically,  Equation  (65)  can  be  put  into  the  same  form  as  Equa¬ 
tion  (62). 

Equations  (55,56,62)  show  that  the  infinitesimal  modulation  response  of  a  soli¬ 
tary  laser  is  externally  manipulated  by  adjusting  the  dimensionless  excess  electrical 
pumping  rate.  Using  the  parameters  in  Table  1,  Figure  4  shows  the  infinitesimal 
modulation  response  of  a  solitary  laser  at  various  steady  state  operating  currents. 
The  function  has  a  Lorentzian  shape  with  a  peak  around  the  relaxation  oscillation 
frequency.  The  trend  in  the  figures  is  that  increased  electrical  pumping  rates  increase 
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the  relaxation  oscillation  and  bandwidth,  but  do  little  to  flatten  the  modulation  re¬ 
sponse. 

The  methods  outlined  in  this  section  for  the  solitary  laser  will  now  be  used  in  the 
model  for  twin  optically  coupled  diode  lasers  to  develop  the  computational  methods  in 
Chapter  III.  The  same  linear  stability  analysis,  outlined  in  this  section,  will  be  applied 
in  Section  3.1  to  determine  the  stability  of  the  steady  states  for  the  coupled  system. 
This  will  lead  into  Section  3.2  where  the  method  for  computing  the  infinitesimal 
modulation  response  of  the  coupled  system  is  derived.  This  computational  method 
will  mirror  the  matrix  treatment  of  the  solitary  laser.  Finally,  a  numerical  method 
will  be  employed  in  Section  3.3  to  calculate  the  nonlinear  modulation  response. 


Figure  4.  Modulation  response  of  a  solitary  InGaAsP  diode  laser  at  an  operating 
current  of  (a)  31.6  mA,  (b)  39.5  mA,  and  (c)  47.4  m A.  The  bandwidths  at  +6  dB  are 
2.4  GHz,  2.9  GHz,  and  3.4  GHz  respectively.  Increasing  current  increases  bandwidth 
but  does  little  to  suppress  the  resonance  at  the  relaxation  oscillation  frequency. 
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III.  Computational  Methods 

In  order  to  compute  the  modulation  response  of  the  coupled  system,  three  computa¬ 
tional  methods  are  developed.  Since  the  diode  lasers  are  modulated  about  a  set  of 
steady  states,  the  stability  of  the  steady  states  needs  to  be  determined  through  a  lin¬ 
ear  stability  analysis.  This  is  accomplished  by  looking  at  small  perturbations  about 
the  steady  states.  If  these  perturbations  grow,  then  the  steady  states  are  unstable  and 
not  appropriate  for  modulation.  Using  concepts  from  the  linear  stability  analysis,  the 
infinitesimal  modulation  response  is  calculated.  This  computation  produces  results 
whether  the  steady  states  are  stable  or  not  and  can  be  used  to  explore  the  parameter 
space.  The  final  computational  method  calculates  the  modulation  response  by  nu¬ 
merically  integrating  the  rate  equations  and  taking  Fourier  transforms.  This  method 
is  used  to  examine  the  modulation  response  at  modulation  amplitudes  beyond  the 
infinitesimal  limit,  but  requires  stable  steady  states. 

3.1  Linear  Stability  Analysis 

The  stability  of  the  steady  states  is  determined  by  examining  the  effect  of  small 
perturbations  about  the  steady  states.  This  is  accomplished  through  a  linearization 
of  the  rate  equations  by  substituting 


Ri(i)  =  RS:  1  +  5Ri(t) 

(66) 

R2(i )  =  Rs,  2  +  §R2(i) 

(67) 

o(i)  =  es  +  so(i) 

(68) 

Mi)  =  Ml  +  $Mi) 

(69) 

Mi)  =  Ml  +  $M(i) 

(70) 

into  Equations  (35-39)  while  taking  P\  and  P2  as  constants  and  ignoring  products 
of  delta  terms.  As  with  the  case  of  the  solitary  laser  in  Section  2.3,  the  result  is 
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expressed  using  the  Jacobian  matrix  of  the  system  as  shown  in  the  following: 


5R1 

5Ri 

5R2 

5R2 

50 

=  J 

50 

J7V, 

5N i 

5N2 

5N2 

(71) 


Solutions  to  Equation  (71)  are  exponential  functions, 


5Ri  (t) 

=  C1AeXli 

+ 

Ch2ex^ 

+ 

C1}3ex^ 

+ 

CXAexd 

+ 

C1Aex^ 

5R2(t) 

=  C2ileXli 

+ 

C%  2ex* 

+ 

C2j3ex* 

+ 

C2Aex 

+ 

C2Aex 

50(t) 

=  C3,  ieAl£ 

+ 

C3,2ex* 

+ 

C3,3ex* 

+ 

C3Aex * 

+ 

C3Aex 

5Nfii) 

=  C4,  ieAl£ 

+ 

CA.2ex^ 

+ 

C4,3ex^ 

+ 

C4Aexd 

+ 

C'4,5eA5,: 

5N2(t) 

=  C5lleAl* 

+ 

C5,2ex* 

+ 

CW3ex^ 

+ 

C5Aex*{ 

+ 

C5,5e^. 

where  A*  are  eigenvalues  of  the  Jacobian  matrix  and  Cm>n  are  arbitrary  constants. 
For  the  system  to  be  stable,  the  real  part  of  all  eigenvalues  must  be  less  than  or  equal 
to  zero.  If  the  real  part  of  any  eigenvalue  is  greater  than  zero  then  there  will  be  a 
term  that  will  grow  exponentially,  thus  destabilizing  the  system. 


3.2  Infinitesimal  Modulation  Response 

The  stability  of  the  steady  states  is  determined  by  examining  small  perturba¬ 
tions  about  the  steady  states,  which  can  also  be  used  to  calculate  the  infinitesimal 
modulation  response.  The  infinitesimal  modulation  response  is  calculated  by  includ¬ 
ing  small  perturbations  in  the  excess  electrical  pumping  rate  along  with  perturbations 
in  the  steady  states.  This  is  accomplished  by  adding  another  vector  to  Equation(71) 
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to  arrive  at 


5R\ 

SR  i 

0 

sr2 

sr2 

0 

SO 

=  J 

SO 

+  SP 

0 

SNi 

SN\ 

1 

T 

sn2 

sn2 

L 

T 

(73) 


where  SP  is  an  infinitesimal  modulation  about  the  steady  state  excess  electrical  pump¬ 
ing  rate.  For  simplicity  and  practicality,  the  lasers  are  modulated  with  the  same  signal 
(SP)  except  for  the  ability  to  modulate  the  signal  out-of-phase  which  is  controlled  by 
the  term  L  [1;  6].  This  term  takes  the  value  of  1  for  in-phase  modulation  and  —1  for 
out-of-phase  modulation.  The  next  step  is  to  take  the  Fourier  transform  of  Equation 
(73)  and  collect  terms  to  arrive  at 


- F{SP } 


0 

FiSR^} 

0 

F{SR2} 

0 

—  [J  iumT\ 

JF{(5©} 

1 

T 

F  {  SN\ } 

L 

T 

F{SN2} 

(74) 


Now,  Equation 


(74)  is  divided  by  F{SP}  and  multiplied  by 


(J-iuI)-1 


to  produce 


F{5R\} 

0 

F{5P} 

FRRz} 

HSP} 

0 

p{se} 

P{SP} 

\J 

0 

HsF} 

1 

P{SP} 

T 

P{SN2} 

L 

L  P{5P}  J 

|_  T  J 

(75) 
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From  Equation  (75)  the  infinitesimal  modulation  response  for  each  laser  is  extracted. 


(76) 


(77) 


where 


H  =  [ J-iuI ]~1. 


The  Matlab®  code  that  implements  this  algorithm  is  included  in  Appendix  C. 


As  long  as  a  set  of  steady  states  can  be  found,  this  algorithm  can  be  used  to  find  the 
infinitesimal  modulation  response  of  a  twin  diode  laser  system  whether  the  steady 
states  are  stable  or  unstable.  The  infinitesimal  modulation  response  provides  some 
information  about  the  response  of  the  system,  but  does  not  provide  information  about 
how  the  system  behaves  when  the  modulation  is  not  infinitesimal. 

A  numerical  algorithm  is  developed  in  the  next  section  to  calculate  the  non¬ 
linear  modulation  response.  This  will  be  used  to  confirm  the  results  of  the  infinites¬ 
imal  modulation  response  as  well  as  provide  a  tool  that  can  be  used  to  examine  the 
modulation  response  when  the  modulation  amplitude  of  the  electrical  pumping  rate 
is  not  infinitesimal. 

3.3  Non-linear  Modulation  Response 

The  modulation  response  is  calculated  numerically  by  first  integrating  the  rate 
equations  for  a  period  of  time  with  a  sinusoidal  modulation  added  to  the  steady  state 
dimensionless  excess  electrical  pumping  rate.  With  this  data,  the  discrete  Fourier 
transform  of  the  slowly  varying  electric  field  amplitude  and  the  discrete  Fourier  trans¬ 
form  of  the  excess  electrical  pumping  rate  are  calculated,  but  only  at  the  frequency 
of  the  sinusoidal  modulation.  This  is  because  the  Fourier  transform  of  the  electrical 
pumping  rate  will  be  a  single  spike  at  the  modulation  frequency.  Essentially,  the 
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sinusoidal  modulation  is  sampling  the  modulation  response  at  one  modulation  fre¬ 
quency.  After  the  Fourier  transforms,  The  absolute  value  of  each  Fourier  transform 
element  is  taken  and  the  appropriate  ratio  is  made.  This  process  is  then  repeated 
for  other  modulation  frequencies  as  desired.  The  C++  source  code  for  the  program 
which  implements  this  algorithm  is  included  in  Appendix  D. 

The  program  starts  by  reading  an  input  hie  which  includes  information  such 
as  laser  parameters,  coupling  parameters,  modulation  frequency  points,  stable  steady 
states  and  variables  to  control  simulation  times.  After  reading  the  input  hie,  the 
integration  is  completed  using  a  fourth  order  Runge-Kutta  algorithm  [24]  with  initial 
values  at  the  steady  states  and  with  a  sinusoidal  modulation  added  to  the  steady 
state  excess  electrical  pumping  rate.  The  Runge-Kutta  algorithm  is 
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and 


£  dRi  .  £  c IR2  .  £  dQ .  r  dNi  .  £  d./V2 

Jl  —  di  >  —  di  5  _  dt  ’  _  di  5  -  dt  * 

To  make  the  discrete  Fourier  transform  element  easier  to  calculate  and  to  ensure 
that  the  modulation  frequency  appears  at  one  of  the  discrete  frequency  points,  the 
integration  time  is  set  to  a  specified  number  of  modulation  periods.  This  also  ensures 
that  the  Fourier  element  of  interest  is  at  the  same  number  of  frequency  points  from 
zero  no  matter  what  the  value  of  the  modulation  frequency.  Once  integration  is 
complete,  the  discrete  Fourier  transform  element  is  calculated  for  the  slowly  varying 
electric  field  amplitudes  out  of  each  laser  diode  and  for  the  excess  electrical  pumping 
rate.  The  Fourier  element  is  calculated  using 

N  N 

F{p  +  !)  =  (0  cos(27rP(^  -  iV-W)  -  i  f(l )  sin(27 rp(l  -  1  )/N),  (79) 

1=1  1=1 

where  F{jp-\- 1)  is  the  Fourier  component,  N  is  the  total  number  of  simulation  points,  / 
is  the  array  that  contains  the  simulation  data  in  the  time  domain,  l  is  the  array  index, 
and  p  is  the  number  of  modulation  periods.  Once  the  Fourier  elements  are  calculated, 
the  absolute  values  of  electric  field  amplitude  Fourier  elements  are  divided  by  the 
absolute  value  of  the  electrical  pumping  rate  Fourier  element.  Finally,  the  results  are 
converted  to  decibels  and  printed  to  the  output  hie.  This  algorithm  calculates  the 
non-normalized  value  of  the  modulation  response  at  one  frequency,  so  the  process  is 
repeated  for  additional  modulation  frequencies.  Outside  of  the  program,  the  results 
are  normalized  to  the  value  of  the  infinitesimal  modulation  response  at  zero  frequency. 
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IV.  Simulations 


In  this  chapter  two  coupling  schemes  are  investigated.  The  first  is  mutual  coupling, 
shown  in  Figure  2(a),  where  photons  are  exchanged  through  the  direct  injection  of 
photons  through  the  output  couplers  [8;  17].  The  second  is  evanescent  coupling, 
shown  in  Figure  2(b),  where  photons  are  exchanged  through  the  evanescent  waves 
that  penetrate  the  waveguides  created  by  the  active  regions  [1;  6;  7;  9;  18;  19]. 

4-1  Mutual  Coupling 

This  section  investigates  injection  coupling  as  shown  in  Figure  2(b)  where  the 
insulating  layer  may  be  a  very  small  air  gap.  For  this  system  a  coupling  phase  of  zero 
is  used  in  Equations  (35-39),  giving 


. jJD 

— =  Aq  R\  +  r]R2  cos(0)  (80) 

dt 

=  N2R2  +  r)R,\  cos(0)  (81) 

dt 

^  =  —a(N2  -  TV,)  -  i j sin(0)  (82) 

TdNj_  =  ^  ^  -  (i  +  2N1)R\  (83) 

dt 

Td-^1  =  p2  -  N2  -  (1  +  2 N2)R22.  (84) 

dt 


Depending  on  the  insulating  material  the  coupling  phase  may  not  be  zero,  but  a  cou¬ 
pling  phase  of  zero  was  chosen  for  ease  of  computing  the  steady  states.  For  simplicity 
the  steady  states  for  Equations  (80-84)  are  calculated  with  symmetric  steady  state 
excess  electrical  pumping  rates  (Ps,i  =  Rs.2)-  It  is  shown  in  Appendix  A  that  the  only 
stable  steady  states  are 

Os  =  0,  Ns  =  -77,  Rs  =  (85) 
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where 


Ns  —  NS)i  —  NS}  2,  Rs  —  Rs,  i  —  Rs,2 


Figure  5  shows  regions  of  stability  for  a  typical  parameter  space. 
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Figure  5.  Regions  of  stability  for  symmetric  steady  states  of  mutual  coupling  with 
symmetric  pumping  rates  and  a  steady  state  phase  difference  of  zero. 


For  a  system  of  twin  InGaAsP  lasers  at  a  normalized  coupling  strength  of  0.1 
and  a  steady  state  current  of  47.4  mA  ( Ps  ps  1),  the  modulation  response  shown  in 
Figure  6  has  a  +6  dB  bandwidth  of  about  3.5  GHz  for  in-phase  modulation  and  a 
—6  dB  bandwidth  of  about  4.4  GHz  for  out-of-phase  modulation.  Figure  6(a)  shows 
that  increased  coupling  for  in-phase  modulation  slightly  increased  bandwidth  as  well 
as  the  position  of  the  peak  response.  Figure  6(b)  shows  that  increased  coupling  for 
out-of-phase  modulation  increases  bandwidth  by  suppressing  the  resonance  at  the 
frequency  of  relaxation  oscillations. 
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(a)  In-phase  Modulation 


Figure  6.  (a)  Modulation  response  of  twin  mutually  coupled  diode  lasers  with 

in-phase  modulation  (b)  Modulation  response  of  twin  mutually  coupled  diode  lasers 
with  out-of-phase  modulation.  The  steady  state  current  is  held  at  47.4  mA.  Solid 
lines  are  the  infinitesimal  modulation  limit  and  dots  are  numerical  results  attained 
by  modulating  the  current  at  an  amplitude  of  316  /j,A . 


^.2  Evanescent  Coupling 


A  system  of  twin  evanescently  coupled  diode  lasers  as  shown  in  Figure  2b  is 
modeled  by  using  a  coupling  phase  of  |  in  Equations  (35-39)  [1;  6]. 

JD 

— J-  =  NiRi  —  rjR.2  sin(0)  (86) 

dt 
rl  R 

=  N2R2  +  gRi  sin(0)  (87) 

dt 

^  =  -a(N2  -Nx)  +  ??cos(0)  (^-  -  (88) 

T<£h  =  _  ^  _  ( 1  +  2Nl)R\  (89) 

dt 

TdN2  =p2_fr2_(1  +  2  n2)R22.  (90) 

dt 

Modeling  the  laser  system  using  Equations  (86-90)  is  discussed  in  References  [1;  6;  18]. 
It  is  shown  in  Reference  [1]  and  Reference  [6]  that  the  laser  system  can  possibly  be 
effectively  modulated  beyond  the  relaxation  oscillation  frequency  by  modulating  the 
lasers  out-of-phase.  A  figure  similar  to  Figure  7a  is  presented  in  Reference  [1]  and 
[6]  showing  the  effects  of  modulating  the  diode  lasers  out-of-phase.  The  eigenvalues 
of  the  Jacobian  indicate  that  there  are  two  resonances.  One  resonance  is  around  the 
frequency  of  relaxation  oscillations  at  about  3.9  GHz  and  the  other  is  at  about  20.5 
GHz.  The  figure  shows  that  the  bandwidth  is  improved  by  about  five  times,  but  the 
steady  states  about  which  the  system  is  modulated  are  unstable.  This  means  that 
the  response  may  never  be  observed.  Appendix  B  shows  the  solution  for  the  steady 
states  when  the  steady  state  excess  electrical  pumping  rates  are  equated  (Ps,i  =  Ps, 2)- 


The  only  real- valued  stable  steady 

states  are 

0s  =  0 

Ns  =  0  Rs  =  \ 

JPs 

(91) 

0S  =  7T 

Ns  =  0  Rs  = 

JPs, 

(92) 
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where 


Ns  —  Ns,  1  —  NS}  2,  —  Rs,l  —  Ra,2, 

but  as  previously  indicated  they  are  unstable  over  a  large  region  of  the  typical  param¬ 
eter  space  as  shown  in  Figure  8.  It  appears  that  the  states  are  stable  when  coupling 
strength  is  very  low  or  when  coupling  strength  is  very  high  and  the  lasers  are  held 
close  to  threshold.  In  the  first  case,  the  lasers  are  essentially  operating  as  solitary 
lasers.  For  the  second  case,  running  the  lasers  close  to  threshold  is  impractical  and 
does  not  fit  well  with  the  model.  This  instability  problem  is  remedied  in  the  next 
sections  by  introducing  asymmetric  pumping  rates  or  gain  saturation  into  the  model. 


Time  (ns) 


(a) 


(b) 


Figure  7.  (a)  Modulation  response  of  twin  evanescently  coupled  InGaAsP  diode 

lasers.  The  steady  state  phase  difference  is  n,  current  through  the  lasers  is  47.4  rnA, 
and  dimensionless  coupling  strength  is  0.1.  The  steady  states  are  unstable  with  the 
secondary  resonance  growing  exponentially  at  a  rate  of  approximately  10los_1.  (b) 
The  time  evolution  of  the  dimensionless  slowly  varying  electric  field  in  one  of  the  lasers 
with  unstable  steady  states  that  are  computed  to  within  10~6  and  no  modulation.  The 
electric  held  starts  at  the  steady  state  value,  but  begins  to  oscillate  at  the  secondary 
resonance  frequency  of  20.5  GHz  and  then  begins  to  oscillate  chaotically. 
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(a)  ©s  =  7 r  Ns  =  0  Rs  =  \J~PS 
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(b)  ©s  =  0  Ns  =  0  Rs=  \/Ps 
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Figure  8.  Regions  of  stability  for  symmetric  steady  states  of  evanescent  coupling 
with  symmetric  pumping  rates  and  a  steady  state  phase  difference  of  (a)  n  and  (b) 
0.  Regions  in  gray  are  where  the  steady  states  are  unstable.  Regions  in  white  are 
regions  where  the  steady  states  are  stable. 
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4-2.1  Asymmetric  Pumping  Stabilization.  For  a  typical  InGaAsP  laser  the 
symmetric  steady  states  with  symmetric  excess  electrical  pumping  rates  are  unstable 
for  a  coupling  strength  of  0.1  and  a  synchronous  steady  state  electrical  current  of  47.4 
mA  (Figure  7).  One  way  of  remedying  this  is  by  making  the  DC  currents  through 
the  lasers  asymmetric.  Asymmetric  steady  states  with  asymmetric  excess  electrical 
pumping  rates  are  found  numerically  using  Newton’s  method  [24],  For  this  algorithm 
the  steady  state  excess  electrical  pumping  rates  are  defined  as 


Ps,2  —  Ps  +  A 


Ps,l  —  Ps  —  A, 


(93) 

(94) 


where  A  is  defined  as  the  dimensionless  excess  electrical  pumping  rate  displacement. 
Substituting  Equations  (93,  94)  into  Equations  (35-39)  yields  a  new  set  of  rate  equa¬ 
tions. 


dU 

— =  N\R\  +  rj R‘2  cos(©  +  0)  (95) 

dt 

dR 

— ir-  =  N2R2  +  r/Ri  cos(©  —  0)  (96) 

dt 

_  =  _a(iv2  -  N,)  -  n  sin(©  -  0)  +  -^  sin(©  +  0)  J  (97) 

pdNi  =ps_A_fri_(1  +  2NJR*  (98) 

dt 

TdN2  =  ps  +  A_jp2_(1  +  2N2)rI  (99) 

dt 

The  term  Ps  is  used  to  calculate  the  symmetric  steady  states  in  Equations  (91,92), 
which  are  then  used  as  the  initial  guess  in  a  Newton’s  method  algorithm  [24],  A  first 
order  correction  to  the  initial  guess  is  made  by  calculating  the  inverse  of  the  Jacobian 
matrix  and  multiplying  it  by  the  rate  equations  evaluated  at  the  initial  conditions. 
This  correction  is  then  subtracted  from  the  initial  guess.  This  process  is  repeated  with 
the  new  steady  state  values  until  the  steady  states  are  reached,  the  maximum  number 
of  iterations  is  exceeded  or  a  tolerance  (Equation  (101))  is  met.  This  algorithm  is 
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shown  in  Equation  (100),  and  the  Matlab®  function  which  implements  it  is  included 
in  Appendix  E. 


Xn-\-l  %r 


j„  f(xr. 


(100) 


where 
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dR2 
dt 
dQ 
dt 
dN i 
dt 

dN2 

dt 


and 


Tol  >  Jn  \f(xn)  •  Jn  \f(xn).  (101) 

By  increasing  the  steady  state  electrical  pumping  rate  displacement  (A),  a  set  of 
stable  steady  states  is  attained  as  shown  in  Figure  9(a).  Figure  10  is  created  using  the 
same  parameters  as  used  to  create  Figure  7,  but  with  a  current  displacement  of  28.4 
rnA.  Since  the  system  is  no  longer  symmetric,  the  modulation  response  out  of  each 
laser  is  different.  Figure  10  shows  that  both  in-phase  and  out-of-phase  modulation 
result  in  a  +6  dB  bandwidth  of  about  4  GHz.  Due  to  the  secondary  resonance  the 
out-of-phase  response  has  an  additional  band  that  spans  about  19.4  GHz  from  about 
4.7  GHz  to  about  24.1  GHz.  This  is  a  significant  improvement  in  bandwidth,  but 
the  different  responses  out  of  each  laser  may  not  be  desirable  in  a  practical  twin  laser 
diode  system.  Stabilization  through  gain  saturation  may  be  a  more  practical  avenue 
for  stability  and  is  addressed  in  the  next  section. 
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12 


(a) 


(b)  (c) 


Figure  9.  Damping  coefficients  vs.  current  displacement.  The  steady  states  become 
stable  at  about  25  mA  of  current  displacement  where  the  coefficient  in  (a)  becomes 
negative. 
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(a)  In-phase  Modulation 


(b)  Out-of-phase  Modulation 


Figure  10.  Modulation  response  for  (a)  in-phase  and  (b)  out-of-phase  modulation  of 
a  typical  InGaAsP  diode  laser.  Solid  lines  are  the  infinitesimal  modulation  limit  and 
dots  are  numerical  results  attained  by  modulating  at  an  amplitude  of  316  /xA.  The 
current  displacement  is  28.4  mA,  the  dimensionless  coupling  strength  is  0.1,  and  the 
phase  difference  is  i r.  The  response  for  in-phase  modulation  has  a  +6  dB  bandwidth 
of  about  4.3  GHz.  The  response  for  out-of-phase  modulation  has  a  +6  dB  bandwidth 
of  about  4.4  GHz,  but  there  is  another  band  that  spans  19.4  GHz  from  the  relaxation 
oscillation  frequency  of  4.7  GHz  to  the  secondary  resonance  frequency  of  24.1  GHz. 
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4-2.2  Gain  Saturation  Stabilization.  Nichols  and  Winful  [19]  show  that  gain 
saturation  stabilizes  the  symmetric  steady  states  shown  in  Equation  (92).  Figure  11 
illustrates  this  finding.  In  order  to  take  advantage  of  this  effect,  a  phenomenological 
two-level  saturable  gain  model  is  used  [19]. 


G(N,E) 


l/rp  +  GN(N-Nth ) 
1  +  /3\E\2 


(102) 


where  0  is  the  gain  saturation  coefficient.  Equation  (102)  replaces  the  linear  gain 
model  in  Equation  (3).  The  rate  equations  are  normalized  using  the  same  variables 
as  shown  in  Equations  (21)  with  the  addition  of  the  dimensionless  gain  saturation 
coefficient  (0). 


P 


w 

GnTc 


(103) 


Typical  values  of  the  gain  saturation  coefficient  are  10~17  cm3  for  InGaAsP  [19]  mak¬ 
ing  typical  values  for  the  dimensionless  gain  saturation  coefficient  around  0.01.  In¬ 
cluding  the  saturable  gain  model,  decomposing,  and  normalizing  yields 


dR\ 

dt 

dR.2 

dt 


I  (rrS  - ^ ^ cos(e  -  ^ 


(104) 

(105) 


dB 

dt 


( 1  +  2  N2 

\1  +  PRI 


1  +  2NX  \ 

i + mi) 


-  V 


sin(0  —  0)  + 


R2 

Ri 


sin(0  +  0) 


(106) 


Tm 

dt 


N\ 


(+±+iW 

\l+PR\) 


(107) 


jdfh 

dt 


P2  -  n2 


1  +  /3  RlJ  - 


(108) 
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The  symmetric  steady  states  with  symmetric  excess  electrical  pumping  rates  are  found 
to  be 


e,=0  r,= 

2  + (3 

O  _  0  PPs  -  2?7  n 

0  s  IT  N g  /S 

2  +  /3 


1  +  2iV,  -  /5(PS  -  i\g 

Ps-Ns 

1  +  2  Ns-p(ps-Ns) 


(109) 

(110) 


By  increasing  the  value  of  /3,  a  set  of  stable  steady  states  is  attained  at  a  value  of 
about  0.02  as  shown  in  Figure  12.  A  gain  saturation  coefficient  of  0.03  was  chosen 
to  calculate  the  modulation  response  shown  in  Figure  13.  The  in-phase  modulation 
response  has  a  —6  dB  bandwidth  of  about  8.2  GHz  and  peaks  around  the  relaxation 
oscillation  frequency  at  3.5  GHz.  The  out-of-phase  modulation  has  a  —6  dB  band¬ 
width  of  about  9.9  GHz  from  the  origin,  but  has  an  additional  —6  dB  band  with  a 
width  of  about  10.5  GHz  from  19.8  GHz  to  about  30.3  GHz.  The  out-of-phase  mod¬ 
ulation  peaks  around  the  secondary  resonance  frequency  of  26.6  GHz.  Comparing 
Figure  7  and  Figure  13  shows  that  a  higher  gain  saturation  coefficient  preserves  the 
basic  features  of  the  modulation  response  when  the  gain  is  linear.  Higher  gain  sat¬ 
uration  seems  to  dampen  the  resonance  at  the  relaxation  oscillation  frequency  while 
enhancing  the  resonance  at  the  secondary  resonance  frequency.  This  stabilization 
scheme  also  does  not  create  two  different  modulation  responses  for  each  laser  as  was 
seen  in  the  asymmetric  stabilization  scheme. 
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Figure  11.  Regions  of  stability  for  symmetric  steady  states  of  evanescent  coupling 
with  symmetric  pumping  rates,  (a)  Steady  state  phase  is  zero  and  gain  saturation  is 
ignored,  (b)  Steady  state  phase  is  tt  and  gain  saturation  is  ignored,  (c)  Steady  state 
phase  is  zero  with  gain  saturation,  (d)  Steady  state  phase  is  tt  with  gain  saturation 
Regions  in  gray  are  where  the  steady  states  are  unstable.  Regions  in  white  are  regions 
where  the  steady  states  are  stable. 
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Figure  12.  Damping  coefficients  vs.  dimensionless  gain  saturation  coefficient.  The 
steady  states  become  stable  at  a  coefficient  value  of  about  0.02  where  the  coefficient 
in  (a)  becomes  negative.  The  split  in  (b)  is  due  to  the  associated  resonance  frequency 
becoming  imaginary. 
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Figure  13.  Modulation  response  of  twin  evanescently  coupled  InGaAsP  diode  lasers. 
Solid  lines  are  the  infinitesimal  modulation  limit  and  dots  are  numerical  results  at¬ 
tained  by  modulating  at  an  amplitude  of  316  /zA.  The  steady  state  phase  difference 
is  7T,  current  through  the  lasers  is  47.4  mA,  dimensionless  coupling  strength  is  0.1, 
and  dimensionless  gain  coefficient  is  0.03.  (a)  The  response  for  in-phase  modulation 
peaks  at  the  free  running  relaxation  oscillation  frequency  of  3.5  GHz  and  has  a  —6  dB 
bandwidth  of  about  8.2  GHz.  (b)  The  response  for  out-of-phase  modulation  peaks  at 
the  secondary  resonance  frequency  of  26.6  GHz  and  has  a  —6  dB  bandwidth  of  about 
9.9  GHz,  but  the  secondary  resonance  creates  another  —6  dB  band  from  about  19.8 
GHz  to  about  30.3  GHz. 
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V.  Summary  and  Conclusions 

Laser  rate  equations  for  electric  field,  and  carrier  density  can  be  used  to  study  the 
modulation  response  of  a  system  of  twin  optically  coupled  diode  lasers.  The  syn¬ 
chronous  steady  states  can  be  calculated  analytically  and  used  in  a  numerical  algo¬ 
rithm  to  calculate  steady  states  where  the  normalized  excess  electrical  pumping  rates 
in  each  laser  are  not  the  same.  The  stability  of  the  steady  states  can  be  examined  us¬ 
ing  a  linear  stability  analysis  that  makes  use  of  the  Jacobian  of  the  rate  equations.  An 
algorithm  using  concepts  from  the  linear  stability  analysis  can  be  used  to  calculate 
the  infinitesimal  modulation  response.  In  addition  to  the  infinitesimal  modulation 
response,  the  non-linear  modulation  response  can  be  calculated  by  using  a  numeri¬ 
cal  method  which  integrates  the  rate  equations,  takes  Fourier  transforms,  and  makes 
the  appropriate  ratios.  For  mutually  coupled  lasers,  it  is  shown  that  the  modula¬ 
tion  bandwidth  of  the  coupled  lasers  does  not  show  significant  improvement  over  the 
modulation  bandwidth  of  the  solitary  laser,  but  the  resonance  at  the  frequency  of  free 
running  relaxation  oscillations  is  suppressed  while  the  coupled  lasers  are  modulated 
out-of-phase  allowing  for  a  flatter  response.  References  [1]  and  [6]  state  that  the  mod¬ 
ulation  bandwidth  for  evanescently  coupled  diode  lasers  is  significantly  larger  than 
the  bandwidth  for  a  solitary  diode  laser  when  modulating  the  signals  out-of-phase, 
but  the  steady  states  about  which  the  lasers  are  modulated  are  unstable.  It  is  shown 
in  this  document  and  Reference  [18]  that  the  steady  states  are  unstable  over  most  of 
the  parameter  space.  This  is  remedied  by  introducing  asymmetric  steady  state  cur¬ 
rents  through  the  laser  or  by  introducing  saturable  gain  [19].  Making  use  of  this,  the 
benefits  of  twin  evanescently  coupled  diode  lasers  as  reported  in  References  [1]  and 
[6]  can  be  realized.  Simulations  in  this  document  show  that  asymmetric  stabilization 
allows  the  coupled  system  to  be  effectively  modulated  out  to  about  25  GHz  with  a  +6 
dB  band  from  0-4.3  GHz  and  a  —6  dB  band  spanning  19.4  GHz  from  4.7-24.1  GHz. 
For  gain  saturation  stabilization,  simulations  show  that  the  coupled  system  can  be 
effectively  modulated  out  to  a  frequency  of  about  25  GHz  with  a  —6  dB  band  from 
0-9.9  GHz  and  another  —6  dB  band  spanning  19.8  GHz  from  10.5-30.3  GHz. 
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There  are  some  additional  analytic,  experimental,  and  modeling  projects  that 
can  grow  form  this  work.  Analytic  expressions  can  be  derived  for  most  of  the  steady 
states  and  infinitesimal  modulation  responses  that  were  calculated  using  the  algo¬ 
rithms  in  Appendix  E  and  Appendix  C,  but  will  require  significant  simplifications 
and  approximations  to  make  them  useful.  Additional  work  can  include  experiments 
were  the  modulation  response  is  calculated  for  evanescent  coupling  [7]  or  for  mutual 
coupling  [8]  and  compared  to  the  results  of  the  model.  Additional  modeling  can 
be  completed  by  investigating  different  coupling  schemes  such  as  changing  coupling 
phases  for  mutual  coupling  and  alternate  configurations  such  as  ring  lasers  [25].  In 
this  work  the  numerical  calculation  of  the  modulation  response  was  used  as  a  con¬ 
firmation  of  the  infinitesimal  modulation  response,  but  can  be  used  to  investigate 
the  non-linear  dynamics  of  the  system.  One  example  of  this  is  the  investigation  of 
sub-harmonic  and  super-harmonic  resonances  [26]  that  occur  when  the  modulation 
amplitude  is  not  infinitesimal  as  shown  in  Figure  14. 


Figure  14.  Non-linear  modulation  response  of  twin  mutually  coupled  InGaAsP 
diode  lasers  with  in-phase  modulation.  Coupling  strength  is  0.001  and  the  steady 
state  current  is  held  at  47.4  m A.  The  legend  indicates  the  modulation  amplitude  of 
the  current. 
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This  work  has  implications  in  current  communications  technology  and  future 
technologies  such  as  chaotic  communications  and  photonic  devices.  Diode  lasers  are 
currently  used  in  todays  communications  systems,  but  there  is  always  room  for  im¬ 
provement.  Increases  in  bandwidth  and  improved  modulation  properties,  realized 
through  coupled  diode  laser  systems,  can  make  current  systems  faster  and  more  ef¬ 
ficient.  It  was  shown  that  the  steady  states  can  be  unstable  for  optically  coupled 
diode  lasers,  but  this  may  not  be  an  unwelcome  circumstance.  These  states  could 
possibly  be  used  to  place  the  system  into  a  chaotic  state  which  can  then  be  used  to 
encrypt  information.  Finally,  this  work  can  be  used  to  drive  upcoming  technologies  in 
optoelectronic  devices  and  photonic  integrated  circuitry,  which  can  lead  to  the  next 
generation  of  communications  and  computing  devices. 
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Appendix  A.  Steady  States  for  Mutual  Coupling 
The  steady  states  of  Equations  (80-84)  where  PSj2  =  Ps,  1  =  Ps  are  solved  by  satisfying 


Ns,iRs,  i  +  yRs,2  cos(0s) 

(Ill) 

Ns,2Rs,2  +  cos(0s) 

(112) 

a(NS}2  Na, i)  77sm(0s)  +  '  ) 

\ris,2  J^s,!  J 

(113) 

Ps  -  Ns,l  -  (1  +  ZNsARh 

(114) 

Ps  ~  Ns, 2  —  (1  +  2Ns,2)R2s  2. 

(115) 

These  equations  are  solved  for  the  steady  states  by  reducing  the  equations  to  a  system 
of  two  equations  with  the  steady  state  dimensionless  excess  carrier  densities  as  the 
variables.  After  this  is  accomplished,  the  new  system  is  solved  for  the  steady  state 
carrier  densities.  The  first  step  is  to  solve  Equation  (111)  for  NSjl,  Equation  (112) 
for  Ns, 2,  Equation  (114)  for  RSt\  and  Equation  (115)  for  Rs>2.  This  result  in 


Ns,  i  = 
Ns,2  = 

Rs,  i  = 

Rs,  2  = 


P'S, 2  { {  \  \ 

-V-fi—  cos(0s) 

ns,i 

Ps,  1  / 

cos(@)s 

-tis,  2 

'Ps  -  Ns,  1 
1  +  2iVs  i 


'Ps-iV. 


s,2 


1  +  2Ns,2 


(116) 

(117) 

(118) 
(119) 


These  equations  show  the  steady  state  carrier  densities  as  functions  of  steady  state 
amplitudes  and  the  steady  state  amplitudes  as  functions  of  steady  state  carrier  den¬ 
sities.  The  phase  dependence  is  eliminated  by  dividing  Equation  (116)  by  Equation 
(117)  to  arrive  at 


Ns,  1  Rj  2 
Ns,  2  Rl  1 


(120) 
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The  first  equation  that  relates  iVSjl  to  iVSi2  is  found  by  substituting  Equations  (118, 
119)  into  Equation  (120)  to  arrive  at 

Ps(Na,2  -  -  (^2  -  till)  -  2 tis,2tiS}1(Ns,2  -  Ns,i)  =  0.  (121) 

In  order  to  find  another  relation  between  NS}i  and  NS}2,  the  trigonometric  identity, 

sin2(@s)  +  cos2(0s)  =  1,  (122) 


is  used.  The  sin2(0s)  is  found  by  solving  either  Equation  (116)  or  Equation  (117) 
for  sin(0s),  squaring  the  result,  and  including  Equation  (120)  where  appropriate. 
The  cos2(@s)  is  found  by  solving  Equation(113)  for  cos(0s),  squaring  the  result,  and 
including  Equation  (120)  where  appropriate.  The  results  of  this  are 


sin2(0s) 


cos2(0s) 


a2(tiSj2  -  Na>iy 


Ns,!NSt2 
r] 2 


(123) 

(124) 


Using  Equation  (122)  and  Equations  (123,124)  a  second  relation  between  iVs>i  and 
1VSj 2  is  found  to  be 


ns,2nsA 


(tis,2  +  iVSii)2  +  a2(lVS)2 


N, 


s,l) 


V2(tis,2  +  Nt 


s,l) 


=  0. 


(125) 


Equations  (121,  125)  are  now  solved  for  the  steady  state  carrier  densities.  For  sim¬ 
plicity 


y  =  Ns,2  +  Ns,i  (126) 

x  =  N,i2  -  NsA  (127) 


45 


are  substituted  into  Equations  (121,125)  to  arrive  at 


0  =  Psx  -  (y2  -  x2)(|)  -  xy 

(128) 

0=  (y  +a  x  )  77  y  . 

(129) 

Equation  (128)  is  now  solved  for  x. 

f 0 

\  ±\J y2  +  2y  -  2Rs- 

(130) 

The  solution  where  x  =  0  is  substituted  into  Equation  (129)  to 

steady  states.  These  steady  states  are 

yield  the  symmetric 

ea  =  0  N.--T,  R,  =  ]j  f*_+^ 

(131) 

@s  =  TT  Ns  —  77  Rs  =  \  — - — , 

(132) 

where 


Ns  —  Ns,i  —  NSj  2,  Rs  —  Rs,  i  —  Rs,  2- 


The  linear  stability  analysis  in  Section  3.1  is  used  to  determine  if  the  steady  states 
are  stable  with  the  parameters  in  Table  1,  a  steady  state  excess  electrical  pumping 
rate  from  zero  to  two  and  a  dimensionless  coupling  strength  from  zero  to  one.  The 
results  of  this  analysis  are  shown  in  Figure  15.  Either  of  the  remaining  solutions  can 
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(a)  (b) 

Figure  15.  Regions  of  stability  for  symmetric  steady  states  of  mutual  coupling 
with  symmetric  pumping  rates  and  a  steady  state  phase  difference  of  (a)  tt  and  (b) 
0.  Regions  in  gray  are  where  the  steady  states  are  unstable.  Regions  in  white  are 
regions  where  the  steady  states  are  stable. 

be  substituted  into  Equation  (129)  to  End 


0  =  y3  (1  +  ct2) 

+  y2  ^2(a2  +  q2)  —  Ps(l  +  a2)^ 

-  y  (4/V) 

+  (2 P,2a2) • 


(133) 


The  asymmetric  steady  states  are  found  by  solving  Equation  (133)  using  the  cu¬ 
bic  formula  found  in  most  algebra  references.  After  solving  for  y,  the  solutions  are 
substituted  into 


Ns,i  =  ^(y  ±  \Jy2  +  2y-  2P^j 

ns, 2  =  \{y^  \Jy2  +  2v- 2 • 


(134) 

(135) 
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to  find  the  steady  state  carrier  densities  which  are  then  used  to  solve  for  the  steady 
state  amplitudes  and  steady  state  phase  difference.  For  these  steady  states  to  be 
valid,  they  must  be  real  and  stable.  For  the  steady  states  to  be  real,  y  must  be  real 
as  well  as  a  part  of  the  space  ((y  <  —1  —  \/l  +  2 Ps)  U  (y  >  — 1  +  \/l  +  2 Ps))  to 
satisfy  Equation  (130).  With  this  and  the  linear  stability  analysis,  the  steady  states 
either  do  not  exist  or  are  unstable  for  the  parameters  in  Table  1,  a  steady  state  excess 
electrical  pumping  rate  from  zero  to  two  and  a  dimensionless  coupling  strength  from 
zero  to  one.  The  results  of  this  analysis  are  shown  in  Figure  16. 
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Figure  16.  Regions  of  stability  and  existence  for  asymmetric  steady  states  of  mutual 
coupling  with  symmetric  pumping  rates.  Regions  in  gray  are  where  the  steady  states 
are  unstable.  Regions  in  dark  gray  are  regions  where  the  steady  states  are  imaginary. 
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Appendix  B.  Steady  States  for  Evanescent  Coupling 

The  steady  states  of  Equations  (86-90)  where  Ps,2  =  Ps,  1  =  E  are  solved  by  satisfying 


NsaRsa  -  gRs,2  sin(0s) 

(136) 

Ns.2Rs.2  +  gRs,i  sin(0s) 

(137) 

a{Ns, 2  Ns,i) +  ycos(Os)  (  ^  ) 

\ris,2  re.5,i  / 

(138) 

Ps  -  Ns,  1  -  (1  +  2 Ns,i)Rh 

(139) 

Ps  —  Ns, 2  —  (1  +  21VS)2)Pg  2- 

(140) 

These  equations  are  solved  for  the  steady  states  by  reducing  the  equations  to  a  system 
of  two  equations  with  the  steady  state  dimensionless  excess  carrier  densities  as  the 
variables.  After  this  is  accomplished,  the  new  system  is  solved  for  the  steady  state 
carrier  densities.  The  first  step  is  to  solve  Equation  (136)  for  NS}1,  Equation  (137) 
for  Ns, 2,  Equation  (139)  for  Rs,i  and  Equation  (140)  for  Rs,2.  This  results  in 


Na,  i  = 

^S,2  •  /y'-X  \ 

V  o  sm(0s) 
ns,i 

(141) 

Ns, 2  = 

,1  •  f  r\\ 

V  „  sm(0)s 
reS)2 

(142) 

Rs,l  = 

IPs  -  Ns,  1 

(143) 

V  1  +  2A^Sii 

Rs,  2  = 

IPs  -  Ns, 2 

(144) 

V  1  +  21VS)2 

This  shows  the  steady  state  carrier  densities  as  functions  of  steady  state  amplitudes 
and  the  steady  state  amplitudes  as  functions  of  steady  state  carrier  densities.  The 
phase  dependence  is  eliminated  by  dividing  Equation  (141)  by  Equation  (142)  to 
arrive  at 


Na,  i 


N. 


s,  2 


Rl  1 


(145) 
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The  first  equation  that  relates  iVSjl  to  iVSi2  is  found  by  substituting  Equations  (143, 
144)  into  Equation  (145)  to  arrive  at 


Ps(Ns,2  +  NSji)  +  4 PsNSj2NsA  -  (1V22  +  N2^)  -  2 NsfiN8ii(Nai2  +  Na>l)  =  0.  (146) 


In  order  to  find  another  relation  between  iVSjl  and  NSj 2,  the  trigonometric  identity, 


sin2(0s)  +  cos2(@s)  =  1, 


(147) 


is  used.  The  sin2(0s)  is  found  by  solving  either  Equation  (141)  or  Equation  (142) 
for  sin(0s),  squaring  the  result,  and  including  Equation  (145)  where  appropriate. 
The  cos2(@s)  is  found  by  solving  Equation(138)  for  cos(0s),  squaring  the  result,  and 
including  Equation  (145)  where  appropriate.  The  results  of  this  are 


sin2(0s) 


cos2(0s) 


NaAN8t2 

7 f 

a2(Ns,2-NS:if 


7]z 


Na, 


NsA 


Nsjl 
Ns, 2 


2) 


(148) 

(149) 


Using  Equation  (147)  and  Equations  (148,  149)  a  second  relation  is  found. 


Na,2NsA  ((Na>2  +  NsA)2  +  a2(NSt2  -  iVSjl)2)  +  r, 2(7VS)2  +  7VM)2  =  0.  (150) 

Equations  (146,  150)  are  now  solved  for  the  steady  state  carrier  densities.  For  sim¬ 
plicity 


y  =  NS)  2  +  NS)  i  (151) 

x  =  Ns, 2  -  N„tl  (152) 
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are  substituted  into  Equations(146,150)  to  arrive  at 


o  =  Psy  +  {y1  -  x2){Ps  -  »-) 

0  =  rfy2  +  - — T—{y2  +  a2x2) 


(153) 

(154) 


Equation  (153)  is  now  solved  for  x. 


x  =  ±1 


'y(i  +  y)(2Ps  -y) 

1  +  2  P8-y 


(155) 


Either  solution  to  x  from  Equation  (155)  can  be  substituted  into  Equation  (154)  to 
arrive  at 


0  =  y5  (1  +  a2) 

-  ?/4  ( (1  —  a2)  +  3PS(1  +  «2)  +  2r/2>) 

V  /  (156) 

-  y3  (4Aa2-(Ps(l  +  «2)+4r?2)(l  +  2Ps)) 

-  y 2  (2r72(l  +  2Ps)2-2P2a2). 

One  solution  for  Equation  (156)  is  y  —  0  which  yields  the  symmetric  steady  states. 
These  steady  states  are 


©s  =  0 

Ns  =  0 

Rs  =  ' 

Jp9 

(157) 

Os  =  7T 

Ns  =  0 

Rs  = 

\fPs, 

(158) 

where 

Ns  —  NSj  i  =  NS)  2,  Rs  =  Rs,  i  =  Rs,  2- 


The  linear  stability  analysis  in  Section  3.1  is  used  to  determine  if  the  steady  states  are 
stable  with  the  parameters  in  Table  1,  a  steady  state  excess  electrical  pumping  rate 
from  zero  to  two  and  a  dimensionless  coupling  strength  from  ICE5  to  one.  The  results 
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of  this  analysis  are  shown  in  Figure  17. 


The  asymmetric  steady  states  are  found 
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Figure  17.  Regions  of  stability  for  symmetric  steady  states  of  evanescent  coupling 
with  symmetric  pumping  rates  and  a  steady  state  phase  difference  of  (a)  tt  and  (b) 
0.  Regions  in  gray  are  where  the  steady  states  are  unstable.  Regions  in  white  are 
regions  where  the  steady  states  are  stable. 


by  solving  Equation  (156)  using  the  cubic  formula  found  in  most  algebra  references. 
After  solving  for  y,  the  solutions  are  substituted  into 


Ah  i  = 


Ah, 2  = 


1/  ,  y(l  +  y)(2  P,-y) 

WW  l  +  2P.-y 


II  ,  y(l  +  y)(2  P,-y) 

^'1  1  +  2  P.-y 


(159) 

(160) 


to  hnd  the  steady  state  carrier  densities  which  can  then  be  used  to  solve  for  the 
steady  state  amplitudes  and  steady  state  phase  difference.  For  these  steady  states  to 
be  valid  they  must  be  real  and  stable.  For  the  steady  states  to  be  real  y  must  be  real 
as  well  as  a  part  of  the  space  {{y  <  —1)  U  {{y  >  0)  fl  (y  <  2 Ps))  U  (y  >  1  +  2 Ps))  to 
satisfy  Equation  (155).  With  this  and  the  linear  stability  analysis,  the  steady  states 
either  do  not  exist  or  are  unstable  for  the  parameters  in  Table  1,  a  steady  state  excess 
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electrical  pumping  rate  from  zero  to  two  and  a  dimensionless  coupling  strength  from 
1CT5  to  one.  The  results  of  this  analysis  are  shown  in  Figure  18. 
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Figure  18.  Regions  of  stability  and  existence  for  asynchronous  steady  states  of 
evanescent  coupling  with  synchronous  pumping  rates.  Regions  in  gray  are  where  the 
steady  states  are  unstable.  Regions  in  dark  gray  are  regions  where  the  steady  states 
are  imaginary. 
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Appendix  C.  ModulatonResponse.m 

This  function  calculates  the  infinitesimal  modulation  response  using  the  procedure  in 
Section  3.2. 

function  [output,  J,  MO]  =  ModulationResponse(varargin) 


if  (strcmp(varargin-[l]- ,  ’help’)) 


fprintf (1 , 
(w,  h,  f, 

fprintf (1 , 
fprintf (1 , 
fprintf (1 , 
fprintf (1 , 
fprintf (1 , 
fprintf (1 , 
fprintf (1 , 
fprintf (1 , 
fprintf (1 , 

fprintf (1 , 
fprintf (1 , 
fprintf (1 , 
fprintf (1 , 
fprintf (1 , 


’[output,  J,  MO]  =  ModulationResponse . . . 
a,  b,  tp,  tc,  L,  rl,  r2,  t,  nl,  n2)\n\n’); 


’ Input 

\n’ )  ; 

’w: 

Normalized  frequencies\n’ ) ; 

’h: 

Coupling  Strength\n’ ) ; 

’f  : 

Coupling  Phase\n’); 

’a: 

Line-width  enhancement  factor\n’); 

’b: 

Gain  saturation  factor\n’); 

’tp: 

Photon  Lifetime  (ps)\n’); 

’tc : 

Carrier  Lifetime  (ns)\n’); 

’L: 

Pump  phase  term  (1  for  in-phase,  -1  for. . . 

out-of-phase)\n’ ) ; 

’rl: 

Steady  state  amplitude  in  laser  l\n’) 

’  r2 : 

Steady  state  amplitude  in  laser  2\n’) 

’t: 

Steady  state  phase\n’); 

’nl: 

Steady  state  carrier  density  in  laser 

l\n’ ) 

’n2 : 

Steady  state  carrier  density  in  laser 

2\n\n 

fprintf (1 , 
fprintf (1 , 
fprintf (1 , 
fprintf (1 , 
return; 


’ Output : \n ’ )  ; 

’output:  Modulation  response\n’ ) ; 

’J:  Jacobian  matrix\n’); 

’MO:  Modulation  response  at  zero  frequency\n\n\n’ ) ; 


end 


w  =  varargin{l} 
h  =  varargin{2]- 
f  =  varargin{3]- 
a  =  varargin{4} 
b  =  varargin{5} 
tp  =  varargin{6} 
tc  =  varargin{7} 
L  =  varargin{8} 
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rl  =  varargin{9} ; 
r2  =  varargin{10} ; 
t  =  varargin{ll} ; 
nl  =  varargin{12} ; 
n2  =  varargin{13} ; 

T  =  tc*1000/tp; 

s  =  size(w) ; 
s  =  s (2) ; 

output  =  zeros (s,  2); 

J  =  [  (1/2) *( (l+2*nl) / (l+b*rl~2)-l) -b* (l+2*nl) *rl~2/ (l+b*rl~2) ~2, 

h*cos (t+f ) , 
-r2*h*sin(t+f) , 
rl/ (l+b*rl~2) , 
0; 

h*cos(t-f ) , 

(l/2)*((l+2*n2)/(l+b*r2~2)-l)-b*(l+2*n2)*r2~2/(l+b*r2~2)~2, 

-rl*h*sin(t-f ) , 
0, 

r2/ (l+b*r2~2) ; 

-h* (sin(t-f )/r2- (r2*sin(t+f ) )/rl~2) -a*b* (l+2*nl) *rl/ (l+b*rl~2) ~2, 
-h* (sin (t+f ) /rl- (rl*sin(t-f ) ) /r2~2) +a*b* (l+2*n2) *r2/ (l+b*r2~2) ~2 , 

-h* (rl*cos (t-f ) /r2  +  r2*cos(t+f )/rl) , 

a/ (l+b*rl~2) , 
-a/ (l+b*r2~2) ; 

(-2* (l+2*nl) *rl/ (l+b*rl~2)+2*b* (l+2*nl) *rl~3/ (l+b*rl~2) ~2)/T, 

0, 

0, 

(-l-2*rl~2/ (l+b*rl~2))/T, 

0; 

0, 

(-2*(l+2*n2)*r2/(l+b*r2~2)+2*b*(l+2*n2)*r2~3/(l+b*r2~2)~2)/T, 

0, 

0, 

(-l-2*r2~2/(l+b*r2~2))/T] ; 


p  =  [0;  0;  0;  1/T;  L/T] ; 
MO  =  abs (inv( J) * (-p) ) ; 
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for  (n  =  l:s) 


W  =  i*w(n)*eye(5) ; 

M  =  abs (inv( J-W) * (-p) ) ; 

output (n, : )  =  [10*logl0(M(l)/M0(l)) ,  10*logl0 (M(2) /MO (2) )] ; 


end 

return; 
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Appendix  D.  NUMOR  Source  Code 

This  program  calculates  the  non-linear  modulation  response  as  outlined  in  Section 
3.3.  A  sample  input  hie  is  shown  below. 

Photon  Lifetime  (ps)  1.600000 

Carrier  Lifetime  (ns)  2.200000 

Linewidth  Enhancement  Factor  5.000000 
Gain  Saturation  Factor 

Coupling  Strength 
Coupling  Phase 

Pumping  Level 
Pumping  Difference 
Modulation  Amplitude 

High  Frequency  (GHz) 

Low  Frequency  (GHz) 

Frequency  Steps 
Frequency  Periods 
Time  Steps 

Steady  States 
rl  1.042955 
r2  1.042955 
t  0.000000 
nl  -0.010000 
n2  -0.010000 

’’Photon  Lifetime”  is  in  picoseconds  and  ’’Carrier  Density”  is  in  nanoseconds.  The 
’’Linewidth  Enhancement  Factor”  is  dimensionless.  The  ’’Gain  Saturation  Factor”  is 
the  dimensionless  gain  saturation  coefficient.  ’’Coupling  Strength”  is  dimensionless 
as  well  as  the  ’’Coupling  Phase”.  The  ’’Pumping  Level”  is  the  steady  state  dimen¬ 
sionless  excess  electrical  pumping  rate  and  the  ’’Pumping  Difference”  is  the  dimen¬ 
sionless  excess  electrical  pumping  rate  displacement.  The  ’’Modulation  Amplitude” 
is  in  dimensionless  form  as  well.  Frequencies  are  in  gigahertz.  ’’High  Frequency”  is 
the  highest  modulation  frequency  to  calculate  and  ’’Low  Frequency”  is  the  lowest 
frequency  to  calculate  with  ’’Frequency  Steps”  as  the  integer  number  of  steps  in  be- 


0.000000 

0.010000 

0.000000 

1.056000 

0.000000 

0.528 

2.500000 

1.000000 

100 

100 

10000 


tween.  ’’Frequency  Periods”  and  ’’Time  Steps”  control  the  integration  times  in  the 
Runge-Kutta  algorithm  and  must  be  integers. 


#def ine 

#include 

#include 

#include 

#include 

#include 

#include 


_CRT_SECURE_NO_DEPRECATE 
<stdio ,h> 

<math . h> 

<conio ,h> 

<time ,h> 

<string.h> 

<stdlib.h> 


//Rate  Equations. 


double  Rlprime (double  n,  double  beta,  double  r,  double  eta,  double  ri, 
double  theta,  double  phi) 

{ 

return  0 . 5* ( (1+2 . 0*n) / (l+beta*pow(r ,2 . 0) ) -1) *r+eta*ri* 
cos(theta+phi) ; 

} 


double  R2prime (double  n,  double  beta,  double  r,  double  eta,  double  ri, 
double  theta,  double  phi) 

{ 

return  0 . 5* ( (1+2 . 0*n) / (l+beta*pow(r ,2 . 0) ) -1) *r+eta*ri*cos (theta-phi) ; 

> 


double  Tprirae (double  alpha,  double  beta,  double  n2,  double  nl, 
double  eta,  double  theta,  double  phi,  double  rl, 
double  r2) 

{ 

return  - (alpha/2 . 0) * ( ( 1+2 . 0*n2) / ( l+beta*pow (r2 , 2 . 0) ) - ( 1+2 . 0*nl) / 
(l+beta*pow(rl ,2.0))) -eta* ( (rl/r2) *sin (theta-phi) + (r2/rl) * 
sin(theta+phi) ) ; 

} 


double  Nlprime (double  p,  double  delta,  double  amp,  double  f, 
double  n,  double  beta,  double  r,  double  T) 


{ 


double  t , 


} 


return  (p-delta+amp*sin(f *t) -n- ( (1+2 . 0*n) / ( l+beta*pow (r , 2 . 0) ) ) * 
pow(r,2.0))/T; 
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double  N2prime (double  p,  double  delta,  double  amp,  double  f,  double  t, 
double  n,  double  beta,  double  r,  double  T,  double  L) 


{ 


> 


return  (p+delta+L*amp*sin(f *t) -n-( (1+2 . 0*n)/ (l+beta*pow(r ,2.0)))* 
pow(r ,2 . 0) )/T; 


// _ 

int  main() 

{ 


//Allocating  memory  to  variables 


FILE*  : 

input ; 

//Pointer  to  input  stream 

FILE*  < 

Dutput ; 

//Pointer  to  output  stream 

double 

pi  =  acos (-1 . 0) ; 

//Pi 

time_t 

ti; 

//Clock  time 

char 

dump [100] ; 

//Character  array  dump 

double 

t; 

//Normalized  time 

double 

dt ; 

//Time  step 

int 

steps ; 

//Number  of  time  steps 

double 

p; 

//Normalized  electrical  pumping  rate 

double 

amp; 

//Amplitude  of  the  electrical  modulation 

double 

delta; 

//P2  =  p  +  delta,  PI  =  p  -  delta 

double 

L; 

//I  for  in-phase  modulation,  -1  for  out 
of  phase  modulation 

double 

eta; 

//Coupling  strength 

double 

phi; 

//Coupling  phase 

double 

taup; 

//Photon  lifetime 

double 

tauc; 

//Carrier  lifetime 

double 

T; 

//Ratio  of  Photon  lifetime  to  carrier 

lifetime 

double 

alpha; 

//Line-width  enhancement  factor 

double 

beta; 

//Gain  saturation  factor 

double 

f; 

//Modulation  frequency 

double 

hif ; 

//High  modulation  frequency  to  calculate 

double 

lof ; 

//Low  modulation  frequency  to  calculate 

double 

bw; 

//High  modulation  frequency  minus  low 
modulation  frequency 

int  periods; 

//Number  of  modulation  periods  to 
calculate 

int 

f steps ; 

//Number  of  modulation  response  points  to 
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double 

reall ; 

double 

imagl ; 

double 

real2 ; 

double 

imag2 ; 

double 

realp; 

double 

imagp ; 

double 

klrl ; 

double 

klr2 ; 

double 

kit ; 

double 

klnl ; 

double 

kln2 ; 

double 

k2rl; 

double 

k2r2 ; 

double 

k2t ; 

double 

k2nl ; 

double 

k2n2 ; 

double 

k3rl ; 

double 

k3r2 ; 

double 

k3t ; 

double 

k3nl ; 

double 

k3n2 ; 

double 

k4rl ; 

double 

k4r2 ; 

double 

k4t ; 

double 

k4nl ; 

double 

k4n2 ; 

calculate 

//Real  part  of  the  Fourier  component  of 
laserl’s  output 

//Imaginary  part  of  the  Fourier  component 
of  laserl’s  output 

//Real  part  of  the  Fourier  component  of 
laser2’s  output 

//Imaginary  part  of  the  Fourier  component 
of  laser2’s  output 

//Real  part  of  the  Fourier  component  of 
the  pump  modulation 

//Imaginary  part  of  the  Fourier  component 
of  the  pump  modulation 

//Variables  for  the  fourth  order 
Runge-Kutta  routine 


// _ 

//Read  input  file  and  initialize  variables _ 

input  =  fopenCinput.dat",  "r"); 
if  (input  ==  NULL) 

{ 

printf ("ERROR:  Input  file  did  not  open\n"); 
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printf ("\nPress  Any  Key  to  Exit"); 


getchO  ; 

return  0; 

f scanf (input , 

"7.s", 

f scanf (input , 

"7.s", 

f scanf (input , 

"•/.If" 

f scanf (input , 

"7.s", 

f scanf (input , 

"7.s", 

f scanf (input , 

"•/.if" 

f scanf (input , 

"7.s", 

f scanf (input , 

"7.s", 

f scanf (input , 

"*/.lf" 

f scanf (input , 

"7.s", 

f scanf (input , 

"7.s", 

f scanf (input , 

"*/.lf" 

T  =  tauc*1000/taup ; 

f scanf (input , 

"7.s", 

f scanf (input , 

"*/.lf" 

f scanf (input , 

"7.s", 

f scanf (input , 

"*/.lf" 

f scanf (input , 

"7.s", 

f scanf (input , 

"*/.lf" 

f scanf (input , 

"7.s", 

f scanf (input , 

"*/.lf" 

f scanf (input , 

"7.s", 

f scanf (input , 

"•/.if" 

f scanf (input , 

"7.s", 

f scanf (input , 

"7.s", 

f scanf (input , 

"*/.lf" 

f scanf (input , 

"7.s", 

f scanf (input , 

"7.s", 

f scanf (input , 

"*/.lf" 

f scanf (input , 

"7.s", 

f scanf (input , 

"7oi" , 

f scanf (input , 

"7oS"  , 

f scanf (input , 

"7oi" , 

f scanf (input , 

"7.S", 

:  (input , 

:  (input , 

:  (input , 
:  (input , 
:  (input , 


"70s",  Mump) 

"7,s" ,  Mump) 

"7oS\  Mump) 
"7.s",  Mump) 

"7oS" ,  Mump) 
"7oS" ,  Mump) 

"7.s",  Mump) 
"7.s",  Mump) 
"7.s",  Mump) 

"7.s",  Mump) 

"7.s" ,  Mump) 

"7.s",  Mump) 
"7.s",  Mump) 
"7oS\  Mump) 
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fscanf (input ,  festeps) ; 


bw  =  hif  -  lof ; 


//Amplitude  of  laser  1  output 

double*  rl  =  (double*) calloc (steps+1 ,  sizeof (double) ) ; 
//Amplitude  of  laser  2  output 

double*  r2  =  (double*) calloc (steps+1 ,  sizeof (double) ) ; 
//Phase 

double*  theta  =  (double*) calloc (steps+1 ,  sizeof (double) ) ; 
//Carrier  density  in  laser  1 

double*  nl  =  (double*) calloc (steps+1 ,  sizeof (double) ) ; 
//Carrier  density  in  laser  2 

double*  n2  =  (double*) calloc (steps+1 ,  sizeof (double) ) ; 
//Pump  modulation 

double*  dp  =  (double*) calloc (steps+1 ,  sizeof (double) ) ; 


fscanf (input , 
fscanf (input , 
fscanf (input , 
fscanf (input , 
fscanf (input , 
fscanf (input , 
fscanf (input , 
fscanf (input , 
fscanf (input , 
fscanf (input , 
fscanf (input , 


"°/0s",  Mump);  fscanf  (input ,  "°/0s",  Mump); 
"°/„s",  Mump); 

"°/.lf\  &rl  [1] ) ; 

"°/„s",  Mump); 

"°/.lf\  &r2  [1]  )  ; 

"°/„s",  Mump); 

"“/.If",  &theta[l] ) ; 

"°/„s",  Mump); 

"°/.lf",  &nl  [1]  )  ; 

"°/„s",  fedump) ; 

"°/.lf",  &n2  [1] )  ; 


f close (input) ; 


// 


//Double  check 


printf ("Photon  Lifetime: 
printf ("Carrier  Lifetime: 
printf ("Linewidth  Enhancement  Factor: 
printf ("Gain  Saturation  Factor: 


°/0f  (ps)  \n" ,  taup) ; 
°/0f  (ns)\n",  tauc)  ; 
%f \n" ,  alpha); 

%f\n\n",  beta); 


printf ("Coupling  Strength: 
printf ("Coupling  Phase: 


0/0f  \n" ,  eta) ; 
7,f  \n\n" ,  phi); 
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°/of  \n" ,  p) ; 

%f\n",  delta); 
°/of\n\n",  amp); 

If  (GHz)\n" ,  hif) ; 

If  (GHz)\n" ,  lof) ; 
°/0i\n" ,  f  steps) ; 

70i\n" ,  periods)  ; 

7oi\n\n" ,  steps)  ; 


printf ("Are  these  correct?  (Y/N) :  "); 
scant  ("7oS" ,  dump); 
if  (strncmp(dump,  "n" ,  1)  ==  0) 

{ 

printf ("\nPress  Any  Key  to  Exit"); 
getchO  ; 
return  0; 

> 

printf ("\n0ut  of  Phase?  (Y/N):  ") ; 
scanf("7oS",  dump); 
if  (strncmp(dump,  "y",  1)  ==  0) 

{ 

L  =  -1.0; 

} 

else 

{ 

L  =  1.0; 

> 


printf ("Pumping  Level: 
printf ("Pumping  Difference: 
printf ("Modulation  Amplitude: 

printf ("High  Frequency: 
printf ("Low  Frequency: 
printf ("Frequency  Steps: 
printf ("Frequency  Periods: 
printf ("Time  Steps: 


printf ("Steady  States\n"); 


printf ("rl: 
printf ("r2: 
printf ("t: 
printf ("nl: 
printf ("n2: 


7of\n" ,  rl  [1] ) ; 

7.f\n",  r2  [1] ) ; 

7«f\n",  theta[l]); 
%f\n",  nl [1] ) ; 
7of\n\n",  n2  [1] ) ; 


// _ 

//Start  the  clock 
time(&ti) ; 
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// 


//Initializing  output  stream _ 

output  =  f open (" output .dat" ,  "w"); 
if  (output  ==  NULL) 

{ 

printf ("ERROR:  Output  file  did  not  open\n"); 

printf ("\nPress  Any  Key  to  Exit"); 
getchO  ; 
return  0; 

} 

f printf (output ,  "Frequency  (GHz)\tModulation  Response  (dB)\t 
Modulation  Response  (dB)\nM); 


// _ 

//Scanning  modulation  frequencies _ 

for  (int  m  =  0;  m  <=  f steps;  m++) 

{ 

printf  ("\nFrequency  Point  70i/°/„i\n",  m,  f steps) ; 

f  =  (lof  +  m*bw/f steps) *taup*0 . 001 ; 
dt  =  periods/(f *steps) ; 
f  =  2.0*pi*f; 
t  =0.0; 

//Fourth  order  Runge-Kutta  routine _ 

for  (int  n  =  2;  n  <=  steps+1;  n++) 

{ 

klrl  =  dt*Rlprime (nl [n-1]  ,  beta,  rl[n-l],  eta, 

r2[n-l],  theta  [n-1]  ,  phi); 
klr2  =  dt*R2prime (n2 [n-1] ,  beta,  r2[n-l],  eta, 

rl[n-l],  theta  [n-1]  ,  phi); 

kit  =  dt*Tprime (alpha,  beta,  n2[n-l],  nl[n-l],  eta, 
theta  [n-1]  ,  phi,  rl[n-l],  r2[n-l]); 
klnl  =  dt*Nlprime(p,  delta,  amp,  f,  t,  nl[n-l],  beta, 

rl  [n-1]  ,  T); 

kln2  =  dt*N2prime (p ,  delta,  amp,  f,  t,  n2[n-l],  beta, 

r2  [n-1] ,  T,  L) ; 
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k2rl  =  dt*Rlprime(nl [n-1]  +  0.5*klnl,  beta, 

rl[n-l]  +  0.5*klrl,  eta,  r2[n-l]  +  0.5*klr2, 
theta [n-1]  +  0.5*klt,  phi); 
k2r2  =  dt*R2prime (n2 [n-1]  +  0.5*kln2,  beta, 

r2[n-l]  +  0.5*klr2,  eta,  rl[n-l]  +  0.5*klrl, 
theta [n-1]  +  0.5*klt,  phi); 
k2t  =  dt*Tprime (alpha,  beta,  n2[n-l]  +  0.5*kln2, 

nl  [n-1]  +  0.5*klnl,  eta,  theta [n-1]  +  0.5*klt, 
phi,  rl [n-1]  +  0.5*klrl,  r2[n-l]  +  0.5*klr2); 
k2nl  =  dt*Nlprime(p,  delta,  amp,  f,  t  +  dt/2, 

nl[n-l]  +  0.5*klnl,  beta,  rl[n-l]  +  0.5*klrl, 
T); 

k2n2  =  dt*N2prime (p ,  delta,  amp,  f,  t  +  dt/2, 

n2  [n-1]  +  0 . 5*kln2 ,  beta,  r2[n-l]  +  0.5*klr2, 
T,  L); 

k3rl  =  dt*Rlprime (nl [n-1]  +  0.5*k2nl,  beta, 

rl  [n-1]  +  0.5*k2rl,  eta,  r2[n-l]  +  0.5*k2r2, 
theta [n-1]  +  0.5*k2t,  phi); 
k3r2  =  dt*R2prime (n2 [n-1]  +  0.5*k2n2,  beta, 

r2 [n-1]  +  0 . 5*k2r2 ,  eta,  rl[n-l]  +  0.5*k2rl, 
theta [n-1]  +  0.5*k2t,  phi); 
k3t  =  dt*Tprime (alpha,  beta,  n2[n-l]  +  0.5*k2n2, 

nl  [n-1]  +  0.5*k2nl,  eta,  theta [n-1]  +  0.5*k2t, 
phi,  rl  [n-1]  +  0.5*k2rl,  r2[n-l]  +  0.5*k2r2); 
k3nl  =  dt*Nlprime(p,  delta,  amp,  f,  t  +  dt/2, 

nl  [n-1]  +  0.5*k2nl,  beta,  rl[n-l]  +  0.5*k2rl, 
T); 

k3n2  =  dt*N2prime(p,  delta,  amp,  f,  t  +  dt/2, 

n2 [n-1]  +  0 . 5*k2n2 ,  beta,  r2[n-l]  +  0.5*k2r2, 
T,  L); 

k4rl  =  dt*Rlprime (nl [n-1]  +  k3nl,  beta, 

rl[n-l]  +  k3rl,  eta,  r2[n-l]  +  k3r2, 
theta [n-1]  +  k3t,  phi); 
k4r2  =  dt*R2prime (n2 [n-1]  +  k3n2,  beta, 

r2[n-l]  +  k3r2,  eta,  rl[n-l]  +  k3rl, 
theta  [n-1]  +  k3t,  phi); 
k4t  =  dt*Tprime (alpha,  beta,  n2[n-l]  +  k3n2, 

nl  [n-1]  +  k3nl,  eta,  theta [n-1]  +  k3t,  phi, 
rl [n-1]  +  k3rl ,  r2[n-l]  +  k3r2) ; 
k4nl  =  dt*Nlprime(p,  delta,  amp,  f,  t  +  dt, 
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nl[n-l]  +  k3nl,  beta,  rl[n-l]  +  k3rl,  T) ; 
k4n2  =  dt*N2prime(p,  delta,  amp,  f,  t  +  dt, 

n2[n-l]  +  k3n2,  beta,  r2[n-l]  +  k3r2,  T,  L) ; 


rl[n] 
r2[n] 
theta  [n] 
nl[n] 
n2  [n] 
dp  [n] 


rl  [n-1]  +  (klrl  +  2.0*k2rl  +  2.0*k3rl  + 
k4rl)/6.0; 

r2  [n-1]  +  (klr2  +  2.0*k2r2  +  2.0*k3r2  + 
k4r2)/6.0; 

theta [n-1]  +  (kit  +  2.0*k2t  +  2.0*k3t  + 
k4t)/6.0; 

nl [n-1]  +  (klnl  +  2.0*k2nl  +  2.0*k3nl  + 
k4nl)/6.0; 

n2  [n-1]  +  (kln2  +  2.0*k2n2  +  2.0*k3n2  + 
k4n2)/6.0; 

amp*sin(f*(t  +  dt)); 


t  =  t  +  dt ; 


printf  ("Percent  Complete:  °/,3.0i\r",  n*100/(steps+l) )  ; 

} 


// 


//Calculate  Fourier  component 


reall  =0.0 
imagl  =0.0 
real2  =0.0 
imag2  =0.0 
realp  =0.0 
imagp  =0.0 


for  (int  1=1;  1  <=  steps+1;  1++) 

{ 


reall  =  reall  + 
imagl  =  imagl  + 
real2  =  real2  + 
imag2  =  imag2  + 
realp  =  realp  + 


rl  [1] *cos(2 . 0*pi*periods* (1-1) / 
(steps+1) ) ; 

rl  [1] *sin(2 . 0*pi*periods* (1-1) / 
(steps+1) ) ; 

r2  [1] *cos(2 . 0*pi*periods* (1-1) / 
(steps+1) ) ; 

r2 [1] *sin(2 . 0*pi*periods* (1-1) / 
(steps+1) ) ; 

dp  [1] *cos(2 . 0*pi*periods* (1-1) / 
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> 


(steps+1) ) ; 

iraagp  =  iraagp  +  dp [1] *sin(2 . O*pi*periods* (1-1) / 
(steps+1) ) ; 


// 


//Print  to  output  stream _ 

fprintf  (output ,  "°/ole\t0/0le\t%le\n" ,  lof  +  m*bw/f steps, 

10 . 0*logl0(sqrt (pow(reall ,  2.0)  +  pow(imagl,  2.0)) 

/sqrt(pow(realp,  2.0)  +  pow(imagp,  2.0))), 
10.0*logl0(sqrt(pow(real2,  2.0)  +  pow(imag2,  2.0)) 

/sqrt(pow(realp,  2.0)  +  pow(imagp,  2.0)))); 


// 

> 


// _ 

//Calculate  run  time,  close  the  output  stream,  and  end  the  program. 

printf ("\n\nThis  run  took  °/0.0f  seconds  to  complete\n" , 
dif f time (time (NULL) ,  ti) ) ; 

f close (output) ; 

printf ("\nPress  Any  Key  to  Exit"); 
getchO  ; 

return  0; 


// 

} 


Appendix  E.  Steady  State,  m 

This  function  calculates  steady  states  using  Newton’s  method  as  outlined  in  Section 

4.2.1.  The  initial  guess  is  the  symmetric  steady  states  with  user  specified  initial 

steady  state  phase  difference.  The  Newton’s  method  algorithm  terminates  when  the 

steady  states  are  reached,  when  the  tolerance  is  met  or  when  the  maximum  number 

if  iterations  is  exceeded.  The  stability  of  these  states  is  determined  using  the  linear 

stability  analysis  in  Section  3.1,  and  the  user  is  inform  if  they  are  stable  or  not. 

Finally  the  user  is  prompted  if  he  wants  to  create  an  input  hie  using  the  calculated 

steady  states  for  the  program  that  calculates  the  non-linear  modulation  response. 

function  [output,  evalues]  =  SteadyState(varargin) 
s  =  size(varargin) ; 
s  =  s  (2)  ; 


if  ((s  ==l)&&(strcmp(varargin-[l} ,  'help’))) 

fprintf(l,  ’[output,  evalues]  =  SteadyState . . . 

(h,  f,  p,  d,  a,  b,  tp,  tc,  t)\n\n’); 
fprintf(l,  ’ Input \n’); 

fprintf(l,  ’h:  Coupling  Strength\n’ ) ; 

fprintf(l,  ’f:  Coupling  Phase\n’); 

fprintf(l,  ’p:  Average  pump  level\n’); 

fprintf(l,  ’d:  Pump  dif f erence\n’ ) ; 

fprintf(l,  ’a:  Line-width  enhancement  factor\n’); 

fprintf(l,  ’b:  Gain  saturation  factor\n’); 

fprintf(l,  ’tp:  Photon  Lifetime  (ps)\n’); 

fprintf(l,  ’tc:  Carrier  Lifetime  (ns)\n’); 

fprintf(l,  ’t:  Steady  state  phase  starting  point\n\n’); 

fprintf(l,  ’ Output : \n ’) ; 

fprintf(l,  ’output:  Steady  States\n’); 

fprintf(l,  ’evalues:  Eigenvalues  of  the  Jacobian  matrix\n\n’ ) ; 
fprintf(l,  ’Tags:\n’); 

fprintf(l,  ’Maxlterations :  Set  maximum  number  of  iterations\n’ ) ; 

fprintf(l,  ’Tolerance:  Set  tolerance  level\n\n’); 

return; 

end 


if  ( (s  ==  10)  |  |  (s  ==  12)  |  |  (s  <  9)  I  I  (s  >  13)) 

fprintf(l,  ’Improper  number  of  input  arguments.  \n’); 
return; 

end 
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h  =  varargin{l}; 
f  =  varargin{2}- ; 
p  =  varargin{3}- ; 
d  =  varargin{4} ; 
a  =  varargin{5} ; 
b  =  varargin{6} ; 
tp  =  varargin{7} ; 
tc  =  varargin{8} ; 
t  =  varargin{9} ; 

T  =  tc*1000/tp; 

if  (s  ==  9) 

tol  =  1.0e-6; 
maxit  =  100; 

end 

if  (s  ==  11) 

if  (strcmp(varargin{10} ,  ’Maxlterations ’ ) ) 
maxit  =  varargindl} ; 
tol  =  10e-6; 

elseif  (strcmp(varargin{10} ,  ’Tolerance’)) 
tol  =  varargindl}; 
maxit  =  100; 

else 

fprintf(l,  ’"7„s"  is  an  improper  tag.  \n’,  varargindO})  ; 
return; 

end 

end 

if  (s  ==  13) 

if  (strcmp(varargin{10},  ’Maxlterations’)) 
maxit  =  varargin{ll} ; 
tol  =  varargin{13} ; 

elseif  (strcmp(varargin{10} ,  ’Tolerance’)) 
tol  =  varargin{ll} ; 
maxit  =  varargin{13} ; 

else 

fprintf(l,  ’  "7„s"  is  an  improper  tag.  \n’,  varargindO})  ; 
return; 

end 

end 
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nl  =  (b*p  -  2*h*cos(t  +  f ) ) / (2  +  b) ; 

n2  =  (b*p  -  2*h*cos (t  -  f ) ) / (2  +  b) ; 

rl  =  sqrt((p  -  nl)/(l  +  2*nl  -  b*(p  -  nl))); 

r2  =  sqrt ( (p  -  n2)/(l  +  2*n2  -  b*(p  -  n2))); 

Rlprirae  =  inline (’ (l/2)*((l+2*n)/(l+b*r~2)-l)*r  +  h*ri*cos (t+f ) ’ , . . . 

’n’ ,  ’b’ ,  >t> ,  ;h’ ,  ’ri\  ’t’ ,  ’f’); 

R2prirae  =  inline(’ (l/2)*((l+2*n)/(l+b*r~2)-l)*r  +  h*ri*cos (t-f ) ’ , . . . 

Jn’,  ’b’ ,  ’r’ ,  ’h’ ,  ’ri’ ,  ’t’ ,  ’f’); 

Tprime  =  inline(,-(a/2)*((l+2*n2)/(l+b*r2~2)-(l+2*nl)/(l+b*rl~2)) • . . 
-h*(rl*sin(t-f )/r2  +  r2*sin(t+f )/rl) ’ , . . . 

’ ,  JbJ ,  J n2J ,  'nl’,  ’h’,  ’t’,  ’f’,  ’rl’,  »r2’); 

Nlprirae  =  inline (’p  -  d  -  n  -  ((l+2*n)/(l+b*r~2))*r~2, , . .  . 

’p\  ’d’,  ’n’,  Jb’,  ’r  ’ ) ; 

N2prirae  =  inline (’p  +  d  -  n  -  ((l+2*n)/(l+b*r~2))*r~2’ , . . . 

;p’>  'd',  ’n’ ,  ’b\  ’r’); 


for  (n  =  l:raaxit) 

J  =  [  (l/2)*((l+2*nl)/(l+b*rl“2)-l)-b*(l+2*nl)*rl“2/(l+b*rl“2)“2, 

h*cos(t+f ) , 
-r2*h*sin(t+f ) , 
rl/ (l+b*rl~2) , 
0; 

h*cos(t-f ) , 

(l/2)*((l+2*n2)/(l+b*r2~2)-l)-b*(l+2*n2)*r2~2/(l+b*r2-2)~2, 

-rl*h*sin(t-f ) , 

0, 

r2/ (l+b*r2~2) ; 

-h* (sin (t-f ) /r2- (r2*sin(t+f ) ) /rl~2) -a*b* (l+2*nl) *rl/ (l+b*rl~2) ~2 , 
-h* (sin (t+f ) /rl- (rl*sin(t-f ) ) /r2~2)+a*b* (l+2*n2) *r2/ (l+b*r2~2) ~2 , 

-h*(rl*cos(t-f )/r2  +  r2*cos(t+f )/rl) , 

a/(l+b*rl~2) , 
-a/ (l+b*r2~2) ; 

(-2* (l+2*nl) *rl/ (l+b*rl~2)+2*b*(l+2*nl)*rl~3/ (l+b*rl~2) ~2)/T, 

0, 

0, 

(-l-2*rl~2/(l+b*rl~2))/T, 

0; 

0, 

(-2*(l+2*n2)*r2/(l+b*r2~2)+2*b*(l+2*n2)*r2~3/(l+b*r2~2)-2)/T, 

0, 
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0, 

(-l-2*r2~2/(l+b*r2~2) )/T] ; 


F  =  [Rlprime (nl ,  b,  rl,  h,  r2,  t,  f ) ; 
R2prime(n2,  b,  r2,  h,  rl,  t,  f); 
Tprime(a,  b,  n2,  nl,  h,  t,  f,  rl,  r2) ; 
Nlprime(p,  d,  nl,  b,  rl) ; 

N2prime(p,  d,  n2,  b,  r2)] ; 


dif  =  inv( [1,0,0,0,0;0,1,0,0,0;0,0,1,0,0;0,0,0,T,0;0,0,0,0,T]*J)*F; 


rl 

=  rl 

-  dif (1) ; 

r2 

=  r2 

-  dif (2); 

t 

=  t  - 

■  dif (3); 

nl 

=  nl 

-  dif (4); 

n2 

=  n2 

-  dif (5) ; 

if 

(sqrt (dif *dif ’ )  <  tol) 

fprintf(l,  ’Steady  states  reached  within  tolerance  at 
iteration:  °/„i  \n\n’ ,  n-1)  ; 
output  =  [rl,  r2,  t,  nl,  n2] ; 
break; 


elseif  (n  ==  maxit) 


output  =  [rl,  r2,  t,  nl,  n2] ; 

fprintf(l,  ’Maximum  number  of  iterations  was  reached  before  ... 
tolerance  was  met.  \n\n’); 


end 


end 


evalues  =  eig(J) ; 


if  ((real(evalues(l)) 
(real (evalues (3) ) 
(real (evalues (5) ) 


<=  0)&&(real(evalues(2)) 
<=  0)&&(real(evalues(4)) 
<=  0)) 


<=  0)&&. . . 
<=  0)&&. . . 


fprintfd,  ’The  steady  states  are  stable\n\n’ )  ; 


else 
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fprintf(l,  ’The  steady  states  are  unstable\n\n’ ) ; 


end 

fprintf(l,  ’Steady  States\n’); 
f printf (1 ,  ’%f\n’,  output); 

str  =  input (’ \nCreate  a  NUMOR  input  file  (y/n) :  ’,  ’s’); 

if  (strcrap(str,  ’y’)) 

fid  =  f open (’ input . dat ’ ,  ’w’); 
if  (fid  ==  -1) 

fprintf (1,  ’\nERR0R:  Input  file  failed  to  open\n’) 
return; 

end 


fprintf (fid, 

’Photon  Lifetime  (ps) 

%f\n’ , 

tp) 

fprintf (fid, 

’Carrier  Lifetime  (ns) 

%f\n’ , 

tc) 

fprintf (fid, 

’Linewidth  Enhancement  Factor 

%f\n’, 

a) ; 

fprintf (fid, 

’Gain  Saturation  Factor 

7of\n\n’  , 

b) ; 

fprintf (fid, 

’Coupling  Strength 

7of\n’ , 

h) ; 

fprintf (fid, 

’Coupling  Phase 

7of\n\n’  , 

f); 

fprintf (fid, 

’Pumping  Level 

7of  \n  ’ , 

p) ; 

fprintf (fid, 

’Pumping  Difference 

7»f  \n  ’ , 

d) ; 

fprintf (fid, 

’Modulation  Amplitude 

\n\n ’ ) ; 

fprintf (fid, 

’High  Frequency  (GHz) 

\n’ ) ; 

fprintf (fid, 

’Low  Frequency  (GHz)  \n 

O; 

fprintf (fid, 

’Frequency  Steps 

\n’ ) ; 

fprintf (fid, 

’Frequency  Periods 

\n’ ) ; 

fprintf (fid, 

’Time  Steps 

\n\n ’ ) ; 

f printf (fid,  ’Steady  States\n’); 
fprintf(fid,  ’rl  °/0f\n’ ,  rl) ; 

fprintf(fid,  ’r2  °/0f\n’ ,  r2) ; 

fprintf(fid,  ’t  %f\n’ ,  t) ; 

fprintf(fid,  ’nl  °/0f\n’  ,  nl) ; 

f  printf  (fid,  ’n2  °/0f  ’  ,  n2) ; 

f  closed  id)  ; 

end 

return; 
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